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ABSTRACT 

In this paper, we describe a new hydrodynamics code for ID and 2D astrophysical simulations, 
BETHE-hydro, that uses time-dependent, arbitrary, unstructured grids. The core of the hydrody- 
namics algorithm is an arbitrary Lagrangian-Eulerian (ALE) approach, in which the gradient and 
divergence operators are made compatible using the support-operator method. We present ID and 
2D gravity solvers that are finite differenced using the support-operator technique, and the resulting 
system of linear equations are solved using the tridiagonal method for ID simulations and an iter- 
ative multigrid-preconditioned conjugate-gradient method for 2D simulations. Rotational terms are 
included for 2D calculations using cylindrical coordinates. We document an incompatibility between 
a subcell pressure algorithm to suppress hourglass motions and the subcell remapping algorithm and 
present a modified subcell pressure scheme that avoids this problem. Strengths of this code include 
a straightforward structure, enabling simple inclusion of additional physics packages, the ability to 
use a general equation of state, and most importantly, the ability to solve self-gravitating hydrody- 
namic flows on time-dependent, arbitrary grids. In what follows, we describe in detail the numerical 
techniques employed and, with a large suite of tests, demonstrate that BETHE-hydro finds accurate 
solutions with 2"'^-order convergence. 
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1. INTRODUCTION 

The ability to simulate hydrodynamic flow is key to 
studying most astrophysical objects. Supernova ex- 
plosions, gamma-ray bursts. X-ray bursts, classical no- 
vae, the outbursts of luminous blue variables (LBVs), 
and stellar winds are just a few phenomena for which 
understanding and numerical tools evolve in tandem. 
This is due in part to the physical complexity, multi- 
dimensional character, and instabilities of such dynam- 
ical fluids. Moreover, rotation is frequently a factor in 
the dynamical development of astrophysical transients. 
One-dimensional hydrodynamics is by and large a solved 
problem, but multi-dimensional hydrodynamics is still a 
challenge. In the context of astrophysical theory, this 
is due primarily to the need to address time-dependent 
gravitational potentials, complicated equations of state 
(EOSs), flexible grids, multi-D shock structures, and 
chaotic and turbulent flows. As a result, theorists who 
aim to understand the Universe devote much of their 
time to code development and testing. 

One of the outstanding and complex problems in 
theoretical astrophysics is the mechanism for core- 
collapse supernovae. For more than two decades, 
the preferred mechanism of explo sion has been 
the delayed neutrino mechanism (jBethe &: WilsonI 
Il985( ). One-dimensional (ID) simulations gener- 
ally f ail to produce exp losions (iLiebendorfer et ail 
2001bllal: rRampp fc Jaiikal 120021: iBuras et all 120031 



Thompson et al.H2003l : rLiebendorfer et al.ll2005f ). How- 
ever, 2-dimensional (2D) simulations, and the ac- 
companying aspherical instabilities, suggest that the 
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neutrino-mechanism may indeed be viable ( Herant et al 
1994; Janka fc Mueller 1995 ; Burrows ct al. 1995, 2007 
Buras et alJ I2OO6I : iKitaura et al. 200 61), though th is has 
yet to be proven. In fact, iBurrows et al.l ()2006l ) have 
recently reported an acoustic mechanism, which seems 
to succeed on long timescales when and if other mech- 
anisms fail. These authors identified two primary rea- 
sons why this mechanism might have been missed before: 
1) 2D radiation- hydrodynamic simulations with reason- 
able approximations are still expensive to run, and with 
limited resources simulations are rarely carried to late 
enough times; and 2), a noteworthy feature of the code 
they used, VULCAN/2D, is its Arbitrary Lagrangian- 
Eulerian (ALE) structure. VULCAN/2D incorporates 
non-standard grids that liberate the inner core from the 
Courant and resolution limitations of standard spherical 
grids. In this context. Burrows ct al. (2006) claim that 
simulating all degrees of freedom leads to a new mecha- 
nism in which the gravitational energy in aspherical ac- 
cretion is converted to explosion energy by first exciting 
protoneutron star core g-modes. These modes then ra- 
diate acoustic power and revive the stalled shock into 
explosion. Remarkably, the acoustic mechanism, given 
enough time, leads to successful core-collapse supernovae 
in all pr ogenitors more massive tha n ~9 M© simulated 
to date (|Burrows et al.ll2006l [2007bD . However, this is a 
remarkable claim, and given the implications, must be 
thoroughly investigated. For i nstance, one q uestion to 
ask is: could the results seen bv IBurrows et al.i (2006.) be 
numerical artifacts of VULCAN? 

Therefore, to address this and other issues that sur- 
round the acoustic supernova mechanism, as well as a 
host of other outstanding astrophysical puzzles, we have 
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designed the new hydrodynamic code, BETHE-hydro"^ , 
from the bottom up. BETHE-hydro wiU be coupled 
wi th the 2D mixed-fr a me ra diation transport scheme 
of iHubenv fc Burrows! ()2007D to create the 2D radia- 
tion/hydrodynamics code BETHE, and the merger of 
these two codes wiU be the subject of a future paper. 
In this paper, we present and test the hydrodynamic 
and gravity algorithms of BETHE-hydro. Since BETHE- 
hydro uses arbitrary grids and a general gravity solver, 
we expect it to be a powerful and flexible numerical tool, 
able to configure the grid to suit the computational chal- 
lenge. 

The core of BETHE-hydro is a ID and 2D ALE hy- 
drodynamics solver. First, solutions to the Lagrangian 
hydrodynamic equations are obtained on an arbitrary 
polygonal grid. Then, to avoid tangled grids, the hy- 
drodynamic variables can be remapped to another grid. 
A unique and powerful feature of ALE schemes is their 
flexibility to tailor an arbitrary grid to the computational 
challenge and to alter the grid dynamically during run- 
time. Hence, purely Lagrangian, purely Eulerian, or ar- 
bitrarily moving grids (chosen to optimize numerical per- 
formance and resolution) are possible. These grids can 
be non-orthogonal, non-spherical, and non-Cartesian. 

Some of the earliest 2D hydrodynamic simulations in 
astrophysics were, in fact, performed to address the 
core-collapse proble m. While some em ployed standard 
fixed-grid schemes (jSmarr et al.l 119811 ). other 2D simu- 
lations were calculated with adaptive grids. Although 
many utilized moving grids, for the most part, the dif- 
feren c ing formulat i ons w e re E ulerian (LcBlanc & Wilson 
flOTOt iSvmbalistvl flOSl iMiller et al.l I1993D . On the 
other hand. iLivio et ahl (jl980D did employ an "Euler- 
Lagrange" method involving a Lagrangian hydrodynamic 
solve, followed by a remapping stage. Even though these 
simulations did exploit radially dynamical grids, they 
were restricted to be orthogonal and spherical. 

Many of the current hydrodynamic algorithms 
used in astrophy s ics ar e based upon either ZEUS 
(jStone fc NormanI Il992f l or higher-order Godunov 
methods, in particu lar the Piecewise- Par abolic 
Method (PPM) (iColella fc Woodward! IT981 
!Woodward fc Colella! !l984f ). Both approaches have 
been limited to orthogonal grids. One concern with 
PPM-based codes is that they solve the hydrodynamic 
equations in dimensionally-split fashion. As a result, 
there have been concerns that these algorithms do 
not adequately resolve flows that are oblique to the 
grid orientation. Rectifying this concern, recent un- 
split higher-order Godunov schemes, or approximations 
thereof, have been develo ped (Truclovc et al. 1998; Klem 
IT999I: iGar^ner fc Stone! l2005t IMiniati fc Colellal 12007^ 
Despite this and the employment of adaptive mesh 
refinement (AMR), these codes must use orthogonal 
grids that are often strictly Cartesian. 

In BETHE-hydro, we use finite-difference 
schemes based upon the su pport-operator m ethod 
dShashkov'fc Steinberg! ! 19951] of !Caramana et al! (|1998! ) 
and !Caramana fc Shashkovr(|1998( ). Differencing by this 
technique enables conservation of energy to roundoff 
error in the absence of rotation and time-varying grav- 

^ Hydrodynamic core of BETHE (Basic Explicit/Implicit 
Transport and Hydrodynamics Explosion Code) 



itational potentials. Similarly, momentum is conserved 
accurately, but due to the artificial viscosity scheme 
that we employ, it is conserved to roundoff error for 
hydrodynamic simulations using Cartesian coordinates 
only. Moreover, using the support-operator method and 
borrowing from an adaptation o f the support-ope rator 
technique for elliptic equations (|Morel et al.! !l998l ). we 
have developed a gravity solver for arbitrary grids. 
Unfortunately, by including a general gravity capability, 
we sacrifice strict energy and momentum conservation. 
However, we have performed tests and for most cases the 
results are reasonably accurate. Furthermore, we have 
included rotational terms and we use a modified version 
of the subc ell pressure scheme to mitig ate hourglass 
instabilities (ICar amana fc ShashkovlfT998fl . 

To resolve shocks, we employ an artificial viscos- 
ity method, which is designed to maint ain grid stabil- 
ity as well (jCampbell fc Shashkov! [200ll ) . and there are 
no restrictive assumptions made about the equation of 
state. Higher-order Godunov techniques employ Rie- 
mann solvers to resolve shocks, but frequently need an 
artificial viscosity scheme to eliminate unwanted post- 
shock ringing. Moreover, the inner workings of Riemann 
solvers often make local approximations that the EOS 
has a gamma-law form, which stipulates that as the in- 
ternal energy goes to zero so does the pressure. For equa- 
tions of state appropriate for core-collapse supernovae, 
this artificially imposed z ero-point energy c an pose prob- 
lems for the simulations (jBuras et al.!!2006h . 

Two codes in astrophysics which have already capi- 
talized o n the arbitrary grid formulations of ALE are 
Djehuty (!Bazan et al.l!2003!: [Dearborn et al.ll2005D and 
VULCAN/2D (|Livne!ll993! : !Livne et al.!!2004l ). In both 
cases, the grids employed were designed to be spheri- 
cal in the outer regions, but to transition smoothly to 
a cylindrical grid near the center. These grid geome- 
tries reflected the basic structure of stars, while avoiding 
the cumbersome singularity of a sphe rical grid. With 
this philosophy, !Dearborn et al.! ()2006! ). using Djehuty, 
have studied the helium core flash phase of stellar evo- 
lution in 3D . In the core-collapse context. !Burrows et al.l 
(|200l[2007ai b). using VULCAN/2D, performed 2D ra- 
diation/hydrodynamic and radiation/MHD simulations 
with rotation. 

Several gravity solvers have been employed in astro- 
physics. The most trivial are static or monopole ap- 
proaches. For arbitrary potentials, the most exten- 
sively used are N-body schemes, which fit most natu- 
rally in Sinooth Particle Hydrodynamics (SPH) codes 
(jMonaghanl !l992l ). As a virtue, SPH solves the equa- 
tions in a grid- free context, and while SPH has opened 
the way to three-dimensional (3D) hydrodynamic simula- 
tions in astrophysics, i ncluding core-collapse simulations 
(jFrver fc Warrenl!2002! ). the smoothing kernels currently 
employed pose serious problems for simulating funda- 
mental hydrodynamic instabilities such as the Kelvin- 
Helmholtz instability (jAgertz et al.! !2007'). Another ap- 
proach to the solution of Poisson's equation for grav- 
ity is the use of multipole expansions of the potential 
(|Miiller fc Steinmetd!l995l ). This technique achieves its 
relative speed by calculating simple integrals on a spher- 
ical grid once. Then, the stored integrals are used in 
subsequent timesteps. In BETHE-hydro, we construct 
finite-difference equations for Poisson's equation using 



the support-operator method and find solutions to the 
resulting linear system of equations for the potentials via 
an iterat ive multiRrid-precondi tioned conjugate-gradient 
method (|Ruge fc Stubenlll987D . 

In Sj2l we give a summary of BETHE-hydro and sketch 
the flowchart of the algorithm. The coordinates and 
mesh details are discussed in SJSI We then describe in 
21 the discrete Lagrangian equations, including specifics 
of the 2"''-order time-integration and rotational terms. 
Section O gives a complete description of the ID and 2D 
gravity solvers. Hydrodynamic boundary conditions arc 
discussed in ^ The artificial viscosity algorithm that 
provides shock resolution and grid stability is described 
in [21 In 3S1 we present the subcell pressure scheme that 
suppresses hourglass modes. Remapping is described in 
depth in fJH In flOl we demonstrate the code's strengths 
and limitations with some test problems. Finally, in JTT] 
we summarize the central characteristics and advantages 
of BETHE-hydro. 

2. BETHE-HYDRO: AN ARBITRARY 
LAGRANGIAN-EULERIAN CODE 

In ALE algorithms, the equations of hydrodynamics 
are solved in Lagrangian form. Within this framework, 
equations for the conservation of mass, momentum, and 
energy are: 

^pV-v, (1) 



Caramana et al.l (I1998D IC aramana k Shashkovl dloMl. 
Campbell fc Shashkovl (|2001f) . and lLoubere fc Shashkovl 
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and 



dt 

dv ^ ^ 

^dt " -P^^-^P' 



de 



(2) 



(3) 



p is the mass density (which we refer to simply as "den- 
sity" and is distinct from the energy or momentum densi- 
ties), V is the fluid velocity, <I> is the gravitational poten- 
tial, P is the isotropic pressure, e is the specific internal 
energy, and d/dt = d / dt + v ■ W is the Lagrangian time 
derivative. The equation of state may have the following 
general form: 

P = P(p,£,X,), (4) 

where Xi denotes the mass fraction of species i. There- 
fore, we also solve the conservation equations: 



dXj 
dt 



0. 



(5) 



Completing the set of equations for self-gravitating as- 
trophysical flows is Poisson's equation for gravity: 

V2$ = 47rGp, (6) 

where G is Newton's gravitational constant. 

All ALE methods have the potential to solve eqs. 
using arbitrary, unstructured grids. In this lies the power 
and functionality of ALE methods. The solutions involve 
two steps and are conceptually quite simple: 1) a La- 
grangian solver moves the nodes of the mesh in response 
to the hydrodynamic forces; and 2) to avoid grid tan- 
gling, the nodes are repositioned, and a remapping algo- 
rithm interpolates hydrodynamic quantities from the old 
grid to this new grid. Of course, the challenge is to find 
accurate solutions, while conserving energy and momen- 
tum. In constructing BETHE-hydro to satisfy these re- 
quirements, we use the ALE hydrodynamic techniques of 



To summarize the overall structure of BETHE-hydro, 
we present a schematic flowchart, Fig. [1] and briefly 
describe the key steps. Establishing the structure for 
all subsequent routines, the first step is to construct the 
arbitrary, unstructured grid. This leads to the next step, 
problem initialization. Then, the main loop for timestep 
integration is entered. After a call to the EOS to obtain 
the pressure, solutions to Poisson's equation for gravity 
are calculated using either the 2D or ID algorithms in 
^ Both use the support-operator method to discretize 
eq. and the resulting system of hnear equations is 
solved by a tridiagonal solver in ID and a multigrid pre- 
conditioned conjugate-gradient iterative method in 2D 

After the timestep is calculated f ii4.2p . the Lagrangian 
equations of hydrodynamics are solved on an arbitrary 
grid using the c ompatible hydrodynamics algorithm of 
I Caramana et al.l (|T99 8) (see 21 for further discussion). 
To ensure 2nd-order accuracy in both space and time, we 
employ a predictor-corrector iteration (i j4.2p . in which a 
second call to the EOS and the gravity solver are re- 
quired. Other than the gravity solver block, the La- 
grangian hydrodynamics solver (Q is represented by the 
set of steps beginning with Predictor and ending with 
Corrector (Fig. [1]). 

After finding Lagrangian hydrodynamic solutions at 
each timestep, one could proceed directly to the next 
timestep, maintaining the Lagrangian character through- 
out the simulation. However, in multiple spatial dimen- 
sions, grid tangling will corrupt fiows with significant vor- 
ticity. In these circumstances, we conservatively remap 
the Lagrangian solution to a new arbitrary grid that mit- 
igates tangled grids. The new grid may be the origi- 
nal grid, in which case we are solving the hydrodynamic 
equations in the Eulerian frame, or it may be a new, 
arbitrarily-defined, grid. 

3. COORDINATE SYSTEMS AND MESH 

In BETHE-hydro, eqs. (HjlH]) may be solved with any 
of several coordinate systems and assumed symmetries. 
Included are the standard ID fc 2D Cartesian coordi- 
nate systems. For astrophysical applications, we use ID 
spherical and 2D cylindrical coordinate systems. Irre- 
spective of the coordinate system, we indicate a position 
by X. The components of the Cartesian coordinate sys- 
tem are x and y; the spherical components are r, 9, and 
(f>; and the cylindrical coordinates are r, z, and 4>. 

Distinct from the coordinate system is the grid, or 
mesh, which defines an arbitrary discretization of space. 
It is this unique feature and foundation of ALE tech- 
niques that provides their flexibility. Figure [2] shows 
two example grids: the butterfly mesh on the left and 
a spiderweb mesh on the right. Either may be used in a 
2D Cartesian or 2D cylindrical simulations, although the 
placement of the nodes is neither Cartesian nor cylindri- 
cal. Instead, these meshes have been designed to sim- 
ulate spherically or cylindrically convergent phenomena 
without the limitation of small zones near the center. 

When using cylindrical coordinates in ALE algorithms, 
one may use control-volume differencing (CVD) or area- 
weighted differencing (AWD); we use CVD. Thorough 
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comparis ons of these two differencing schemes may be 
found in ICaramana et aTl ([1998). Here, we give a basic 
justification for choosing CVD. For CVD, the volumes 
of the subcell, cell, and node are calculated by straight- 
forward partitioning of these regions by edges. Hence, 
volumes and masses are exact representations of their re- 
spective regions. While this discretization is natural and 
easy to comprehend, it does not preserve strict spherical 
symmetry when using a spherical grid and cylindrical co- 
ordinates. On the other hand, AWD is designed to pre- 
serve spherically symmetric flows, but only with the use 
of spherical grids that have equal spacing in angle. Since 
this prevents the use of arbitrary grids, and the asymme- 
tries of CVD are small, we use CVD for discretization. 

Construction of the grid begins with the arbitrary 
placement of nodes. At these nodes, coordinate posi- 
tions, velocities, and accelerations are defined. Simply 
specifying node locations does not completely define a 
mesh, since there is not a unique way to assign cells 
and masses to these nodes. Therefore, the user must 
specify the connectivity among nodes, which define the 
arbitrarily-shaped polygonal cells. It is within these cells 
that cell-centered averages of p, e, P, Xi, and $ are de- 
fined. 

With the node positions, their connectivity, and the 
cells defined, a mesh is completely specified and all other 
useful descriptions follow. Cells are denoted by z, and 
nodes are indexed by p. The set of nodes that defines 
cell z are p £ S{z), where the nodes are ordered counter- 
clockwise. Conversely, the set of cells that shares node p 
is denoted by z S S{p). Each cell has Np{z) nodes that 
define it, and each node has Nz (p) cells that share it. The 
sample sub- grid depicted in Fig. [3] helps to illustrate the 
nodal and cell structure. The filled circles indicate node 
positions, Xp, and the crosses indicate the cell-center po- 
sitions, Xz = {'^/Np{z))J2p£S{z) ^p- solid lines are 
direction-oriented edges that separate cells from one an- 
other, and the open circles indicate their mid-edge loca- 
tions. Partitioning the cell into subcells, the dashed lines 
connect the mid-edges with the cell centers. No matter 
how many sides a cell has, with this particular division 
the subcells are always quadrilaterals. Naturally, the cell 
volumes (see appendix|B]for formulae calculating discrete 
volumes), Vz, are related to the subcell volumes, Vp, by 

Furthermore, each node has a volume, Vp, defined by the 
adjoining subcells that share the node p: 

Vp^ J2 (8) 

zesip) 

For calculating pressure forces and fiuxes, vectors are as- 
signed to each half edge on either side of node p: Cp_^^ 

and Cp_, where -I- indicates the half edge in the coun- 
terclockwise direction around the cell and — indicates 
its opposite counterpart. Their magnitudes are the ar- 
eas represented by the half edges, and their directions 
point outward and normal to the surface of cell z. From 
these half-edge area vectors, an area vector, Sp, that is 
associated with zone z and node p is then defined: 

Sp = C^p+ + . (9) 



Similarly, vectors are associated with the lines connect- 
ing the mid-edges and the cell centers: dp. Again, their 
magnitudes are the corresponding area, but while the 
vector is normal to this line, the direction is oriented 
counterclockwise around the cell. 

4. DISCRETE LAGRANGIAN HYDRODYNAMICS 

The fundamental assumption of Lagrangian algorithms 
is that the mass, m^, of a discrete volume Vz is constant 
with time. For staggered-grid methods, in which scalars 
are defined as cell-centered averages and vectors are de- 
fined at the nodes, it is necessary to define a Lagrangian 
mass, nip, for the nodes as well. This nodal mass is 
associated with the node's volume Vp (eq. [8|). Conser- 
vation of mass (eq. [1]) implies zero mass flux across the 
boundaries, dV, of either the cell volume or nodal vol- 
ume. Therefore, the region of overlap for Vz and Vp, 
which is the subcell volume Vp, is bounded by surfaces 
with zero mass flux. Consequently, the most elemental 
Lagrangian mass is the subcell mass, m^, and the cell 
mass (niz) and node mass (mp) are constructed as ap- 
propriate sums of subcell masses: 

mz= ^P' (10) 

and 

mp= J2 (11) 

2es(p) 

Hence, we arrive at three discrete forms of mass conser- 
vation (eq. [T|): 

p;^^, Pp^^, and pz^^. (12) 

In defining the discrete momentum and energy equa- 
tions, we use the compatible h y drody namics algorithms 
developed by ICaramana et all (|l998i ). Specifically, the 
discrete divergence and gradient operators are compati- 
ble in that they faithfully represent their analog in con- 
tinuous space and their definitions are expressly related 
to one another using the hydrodynamic expressions for 
conservation of momentum and energy. As a result, 
this approach leads to discretizations that satisfy mo- 
mentum and energy conservation to machine accuracy. 
This is accomplished with th e support-operator method 
(jShashkov fc Steinbergl [r995f l . Given the integral iden- 
tity, 

/ <i>{W-H)dV+ I H-V<^dV=<h ^H-dS, (13) 
Jv Jv Jav 

where $ is any scalar and H is some vector, there is an 
incontrovertible connection between the divergence and 
gradient operators. For many choices of discretization, 
the discrete counterparts of these operators could vio- 
late this integral identity. Simply put, the goal of the 
support-operator method is to define the discrete oper- 
ators so that they satisfy eq. (fT3|) . The first step is to 
define one of the discrete operators. It is often, but not 
necessary, that the discrete divergence operator is defined 
via Gauss's Law: 

/ ^ -vdV ^ f v-dS, (14) 

Jv JdV 
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and then eq. p3)) is used to compatibly define the other 
discrete operator. Discretizing the hydrodynamic equa- 
tions, ICar^^an^^gOn 



begin by defining (VP)p at 
the nodes and use the integral form of energy conserva- 
tion and an equation equivalent to eq. (fT3|) to define for 

each cell the discr ete divergence of the velocity, (V • v)z- 
To begin, C aramana et al.l (|1998D integrate eq. ^ (ex- 
cluding gravity and rotational terms) over the volume of 

node p, producing the discrete form for (VP)p and the 
momentum equation: 



' At 



VPdV 



PdS 



dV„ 



where An 



p 



V, 



P 



is the change in velocity from 
timestep n to the next timestep n + I, the timestep is 
At, and is the pressure in cell z. In other words, the 
net force on node p is a sum of the pressure times the 
directed zone areas that share node p. Hence, the subcell 
force exerted by zone z on point p is = PzS^. A more 
complete description of the subcell force, however, must 
account for artificial viscosity (fjT]): 



Fz=pJl + ft,, 



(16) 



where /^'visci is the subcell force due to artificial viscos- 
ity. Furthermore, in this work, we include gravity and 
rotation for 2D axisymmetric simulations, and the full 
discrete momentum equation becomes 



Avp 
At 



E /I 

zes(p) 



nipgp 



uipAp , 



(17) 



where gp and Ap are the gravitational and rotational 
accelerations, respect ively. For simpl i city, a nd to par- 
allel the discussion in ICaramana et al.1 ()1998[ ) , we ignore 
these terms in the momentum equation and proceed with 
the compatible construction of the energy equation (see 
i j4.3.1l and ii5.4l for discussions of total energy conserva- 
tion including rotation and gravity, respectively). 

To construct the discrete energy equation, 
ICaramana et al.l ()1998[ ) integrate eq. ([3|) over the 
discrete volume of cell z: 

Ae P ^ 

= - PV • vdV = -P,(V • v)zVz , (18) 

where Asz = £^~^^ — Then, the objective is 

to determine a discrete form for the right hand side 
of this equation that conserves energy and makes the 
discrete g radient and divergence operators compatible. 
ICaraman a et al.l (jl998l ) accomplish this with the integral 
for conservation of energy (neglecting gravity) : 



de 1 dv'^ 



dt 



2^ dt 



dV-^ 



{PV ■v + v-VP)dV 



Pv-dS, 



(19) 



where the second expression is the integral identity, eq. 
p^ . that defines the physical relationship between the 
gradient and divergence operators. Neglecting boundary 
terms, the discrete form of this integral is 



^1 At 



pes{z) 



(20) 



where we used Avp — [v^'^ ) — [v^) — 2vp ' ■ Avp, 
^ i/2{v^+'^+v^), and substituted in eq. If 
we set the expression for each zone in eq. ((20)) to zero, 
we arrive at the compatible energy equation: 

-^^ = - E -T'^'-fl- (21) 

pesiz) 



The RHS of eq. pTjl is merely the PdV work term of the 
first law of thermodynamics. Its unconventional form is 
a consequence of the support-operator method. 

Thus, we derive two significant results. First, inspec- 
tion of eqs. (PT|) and P3)) leads to compatible defini- 
tions for the discrete gradient and divergence operators. 
Specifically, the discrete analogs of the gradient and di- 
vergence operators are 



and 



P zesip) 



(22) 



(23) 



peS{z 



Second, eqs. ^5]) . and (PT|) form the compatible 

discrete equations of Lagrangian hydrodynamics. 

4.1. Momentum Conservation 

In its current form, BETHE-hydro strictly conserves 
momentum for simulations that use Cartesian coordi- 
nates and not cylindrical coordinates. ICaramana et all 
show that the requirement for strict conservation 



is the use of control- volume differencing, or that exact 
representations of cell surfaces are employed. As a simple 
example, consider momentum conservation with pressure 
forces only: 



E^p-E E /1 = E^^ E Sp 



(24) 



zes{p) 



pes{z) 



dV 



For Cartesian coordinates, the sum X)pGS(z) ^p exactly 
zero for each cell as long as the surface-area vectors are 
exact representations of the cell's surface. This is not the 
case for area-weighted differencing, and hence momen- 
tum is not conserved in that case. For cylindrical coordi- 
nates, the sum of surface-area vectors gives zero only for 
the z-component. However, the assumption of symmetry 
about the cylindrical axis ensures momentum conserva- 
tion for the r-component. Therefore, control-volume dif- 
ferencing ensures momentum conservation, even with the 
use of cylindrical coordinates. 

Based upon similar arguments, the subcell-pressure 
forces that eliminate hourglass motions (^5]) are formu- 
lated to conserve momentum for both Cartesian and 
cylindrical coordinates. On the other hand, because we 
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multiply the Cartesian artificial viscosity force by 2nrp to 
obtain the force appropriate for cylindrical coordinates, 
the artificial viscosity scheme (fJ7]) that we employ is con- 
servative for Cartesian coordinates only. 

We test momentum conservation with the Sedov blast 
wave problem (i jlO.Sp using Cartesian and cylindrical 
coordinates. Symmetry dictates conservation in the r- 
component of the momentum. Consequently, we com- 
pare the z-component of the total momentum at t = 
s with subsequent times. Since the initial momentum of 
this test is zero, we obtain a relative error in momentum 
by Pz{t)/Pz,z>o{t), where Pz{t) is the z-component of 
the total momentum at time t and Pz,z>o{t) is a similar 
quantity for z > 0. As expected, the relative error using 
Cartesian coordinates reflects roundoff error. For cylin- 
drical coordinates and 35,000 cells, the relative error in 
momentum is ~ 3 x 10^^ after 0.8 s and 71,631 timesteps. 
Although momentum is strictly conserved for simulations 
using Cartesian coordinates only, the Sedov blast wave 
problem demonstrates that momentum is conserved very 
accurately, even when cylindrical coordinates are used. 

4.2. Second- Order Time Integration: 
Predictor- Corrector 

Our method is explicit in time. Therefore, for numeri- 
cal convergence, we limit the timcstep, Ai, to be smaller 
than three timescales. By the Courant-Friedrichs-Levy 
(CFL) condition we limit At to be smaller than the short- 
est sound-crossing time and the time it takes flow to tra- 
verse a cell. The latter condition is useful to avoid tan- 
gled meshes and is necessary when remapping is used. In 
addition, it is limited to be smaller than the fractional 
change in the volume represented by {S/ ■v)z- In calculat- 
ing these timescales a characteristic length for each cell is 
needed. Given arbitrary polygonal grids, we simply use 
the shortest edge for each cell. In practice, the timestep 
is the shortest of these timescales multiplied by a scaling 
factor, CFL. In our calculations, we set CFL = 0.25. 

To ensure second-order accuracy and consistency with 
the time levels speciflcd in eqs. ([T7| and ((2T|) . we em- 
ploy predictor-correct or time integratio n using the finite 
difference stencil from ICaramana et al.l (fT998>) . For each 
timcstep, the current time level of each quantity is iden- 
tified with n in superscript. The time-centered quantities 
are indicated with a n+ 1/2 superscript, while the pre- 
dicted values are labeled n-\- l,pr, and the fully updated 
values for the next timestep are labeled with a n -t- 1 
superscript. 

Table [1] lists the steps for the predictor-corrector in- 
tegration. The predictor step begins by defining the 
forces on nodes, /", using current pressures P", areas 
5", and artificial viscosity forces f^^^c- ^q. (fTT]) is 
then used to calculate the predicted velocity, v^+'^^p^ , 
The predicted node coordinate is then calculated by 
^+i..pr = ^ + Aiir"+i/2, where v^+'^l^ is the time- 
centered velocity computed as — i^^i+i^P'' +7/"). 

Then, using /" and v^+'^/'^ in eq. (|2T|) . the predicted spe- 
cific internal energy, £"+1'P''^ is obtained. Using the pre- 
dicted node positions, the predicted rotational accelera- 
tion volume (l/«+i^P'^), and density 
are computed, which in turn are used to calculate the 
predicted gravitational acceleration, g^+'^-P^, Calling the 



EOS with e"+i.P'" and p^+^^P"^ gives the predicted pres- 
sure We finish the predictor step and prepare 
for the corrector step by calculating the time-centered 
positions x^'^'^^'^, pressures P"+i/2^ rotational accelera- 
tions ^p"*"^^^, and gravitational accelerations (f^+^Z^ as 
simple averages of their values at the {n -\- l,pf) and n 
time levels. 

The corrector step then proceeds similarly to the pre- 
dictor step, and 5'"+i/2 are used to compute 
which in turn is used to compute the updated 
velocity, v^'^^. After computing the new time-centered 
velocity, the node coordinates are finally updated using 
the time-centered velocity. Then, the specific internal 
energy, e"+^/^, is updated using and Fi- 
nally, the volume, and density are computed 
for the next timestep. 

4.3. Rotation in 2D cylindrical coordinates 

In 2D simulations using cylindrical coordinates, we 
have extended BETHE-hydro to include rotation about 
the axis of symmetry. For azimuthally symmetric con- 
figurations, all partial derivatives involving (jj are zero. 
Consequently, the conservation of mass and energy equa- 
tions are 

-p^rz ■ V (25) 



and 



Dt 

De 



-PVr 



(26) 



(27) 



where the rz subscript on the del operator reminds us 
that the derivatives with respect to 4> vanish. For concise 
notation, we define a pseudo Lagrangian derivative: 

D _ d d d 

While eqs. and are simple extensions of the true 
mass (eq. [1} and energy (eq. [3]) equations, the presence 
of "Christoffel" terms in the momentum equation give it 
extra terms: 

P^^-VrzP + pA. (28) 

where the extra rotational terms are represented by A: 

OBz 

A={ < e.. (29) 



vJv,. 



To see the origin of these extra terms, one must begin 
with the momentum equation in cylindrical coordinates 
and Eulerian form: 

\ at ^ dz ^ Qj. + ^ J - 

+ + + + = -^f • 

(30) 

Again, dropping terms involving partial derivatives of 0, 
and moving all terms not involving partial derivatives to 
the right-hand side (RHS), eq. ([30|) reduces to 



Vz 


dvz 
dz 


f Vr 


dvz 
dr 


Vz 


dVr 

dz 


f Vr 


dvj. 
dr ^ 


Vz- 


dz 


- Vr'- 


dv^ 
dr 



dP 

or r 



(31) 
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Hence, using our definition for the pseudo Lagrangian 
derivative, D/Dt, and V^z in eq. ([?T|) we obtain eq. (pS)) . 
It is now apparent why D/Dt is a pseudo Lagrangian 
derivative. Strictly, the "extra" rotational terms are a re- 
sult of "Christoffel" terms in the true Lagrangian deriva- 
tive. Therefore, in defining D/Dt without them, we've 
defined a pseudo Lagrangian derivative. 
Conveniently, the component of eq. 



Dv^ 

-) 

Dt 



-P 



V^Vr 



(32) 



is simply a statement of the conservation of angular mo- 
mentum about the axis of symmetry: 



DJ 
Ut 



= 0. 



(33) 



where J is the angular momentum about the z-axis of a 
Lagrangian mass. 

Other than the appearance of an extra term on the 
RHS of the r-component and of an equation for the (j) 
component in eq. (|28p . the 2D hydrodynamic equations 
including rotation about the z-axis are similar to the 
equations without rotation. Therefore, the algorithms 
for dynamics in the 2D plane remain the same, except for 
the inclusion of the conservation of angular momentum 
equation and the fictitious force term. 

It seems natural to define J — TripV^^prp and v^ p — 
LOpVp, where tUp is the angular velocity of node p. In 
fact, our first implementation of rotation did just this. 
However, results indicated that this definition presents 
problems near the axis of rotation. Whether u^.p or ujp 
is defined, the angular momentum at node p is associ- 
ated with the region of subccUs that comprise the node's 
volume. For the nodes on the axis, this definition states 
that J = since rp = 0, which implies that the angu- 
lar momentum for the region of subcells near the axis 
is zero. This awkward feature causes no problems for 
Lagrangian calculations with no remapping. However, 
during the subcell remapping step, the subcells near the 
axis contain no angular momentum, yet physically there 
should be some angular momentum to be remapped. 

Therefore, we have devised an alternative method for 
including rotational terms in 2D simulations. Just as 
for mass, the angular momentum of a subcell, J^, is a 



primary indivisible unit, which by eq. (|33p is constant 
during a Lagrangian hydrodynamic timestep. The angu- 
lar momentum for a cell is 



J. 



E 

pes(z) 



J' 



and for each node is 



(34) 



(35) 



To relate angular velocity and to the angular mo- 
mentum, we treat each subcell as a unit with constant 
angular velocity. For such regions of space the angular 
momentum is related to the angular velocity by 



p p 



(36) 



where Ip = Jy^ PpT^dV is the moment of inertia for sub- 
cell {z,p}. Since we are interested in accelerations and 



velocities at the nodes, we choose that ojp 
naturally leads to 



^plp 



which 



(37) 



Substituting eq. (I37p into eq. ([35|) . we define a relation- 
ship between the angular momentum at the nodes and 
the angular velocity at the nodes: 

(38) 
(39) 



where 



zesip) 



By conservation of angular momentum, Jp^^ — Jp = 
Jp — ti>"+^X"+-'^, the angular velocity at the n+1 timestep 



n+l _ Jp 



rn+l 



(40) 



Finally, a new form for the r-component of Ap is con- 
structed. The previously defined form is v'^/r, and hav- 
ing r in the denominator poses numerical problems near 
and at r = 0. Since v'^/r ~ w^r, the new half timestep 
acceleration is 



Ar 



n+l/2 



(41) 



4.3.1. 



Conservation of Energy with Rotation 

Without rotation and time-varying gravity, the equa- 
tions for 2D simulations using cylindrical coordinates 
conserve energy to machine precision. However, with 
rotation, strict energy conservation is lost. Using the 
equations for energy and momentum including rotation, 
eqs. (l26l) and ([28]) the energy integral becomes 



De 1 Dv'^ 



dV = 



-PVr 



-V-VrzP + PiVrAr+V^A^)j dV . (42) 

Since the terms VrAr and Vff,A^ are equal in magnitude 
and opposite in sign, the rate of change for the total 
energy reduces to the usual surface integral: 



dE 
'dt 



— = <h Pv-dS. 



(43) 



If the form of the energy equation is to remain the same, 
discrete versions of VrAj. and v^A^ must be derived that 
cancel one another. 

In the previous section, we made arguments not to use 
J = rripV^r. However, it is simpler and more intuitive to 
illustrate why the 2D equations, including rotation, can 
not satisfy strict energy conservation, while satisfying 
angular momentum conservation. With conservation of 
angular momentum we have ■y^+ir""'""'^ = ti^r" or 



The discrete time derivative of is then 



n n+l/2 
V^Vr 



At 



„ji+i 



(44) 



(45) 
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where we have used Ar = w"^^ Ai. This specifics the 
required time levels for the components of to conserve 
angular momentum. Considering time integration, the 

discrete analog of VrAr + v^A^ is Vr' 

To ensure that this sums to zero, the discrete analog of 

Ar is 



n+l/2 n 



(46) 



With these time levels, the 2D equations, which in- 
clude rotation, would satisfy energy and angular momen- 
tum conservation to round-off error. However, r""*"^ and 

v^^'^^^^ are not known at the beginning of the timestep. It 
is only after the predictor step of the predictor-corrector 
method that v'^'^^^'^ and are obtained. As a re- 

sult, Aj., including the correct time levels and predictor 
values, is 



A-r — 



n+l/2 „ 



(47) 



but to conserve angular momentum and energy, r"^^ 
appears in the denominator. Therefore, the terms 

v^'^^^'^Acf, and Vr^^^'^ Ar will not be exactly opposite and 
equal in magnitude, and angular momentum and energy 
conservation can not be enforced to machine precision at 
the same time. 

In the previous section, we motivated a different dis- 
cretization of angular momentum, angular velocity, and 
However, the basic results of this section do not 
change. There is still the fundamental problem that 
to strictly conserve angular momentum, positions at the 
n + 1 time level are required, while one must use the pre- 
dicted, n + l,pr, values in updating the radial velocity. 
Sacrificing strict energy conservation, we have chosen to 
strictly conserve angular momentum for simulations in- 
cluding rotation. 

4.3.2. Simple Rotational Test 

We have designed a test that isolates and tests rota- 
tion. Using cylindrical coordinates, the calculational do- 
main lies in the r — z plane and extends from the sym- 
metry axis, r = 0, to r = 1.0 cm and from z = —0.1 
cm to z = 0.1 cm . The initial density is homogeneous 
with p = 1.0 g cm~^. Isolating the rotational forces, we 
eliminate pressure gradients by setting the internal en- 
ergy everywhere to zero. Consequently, the trajectory 
of a Lagrangian mass is a straight line in 3D and has a 
very simple analytic evolution for r(t), determined only 
by its initial position and velocity. We set — and 
give U0 and Vr homologous profiles: Vr = —3r/ro cm 
and U0 — r/ro cm s~^. This produces a homologous 
self-similar solution for Vr and that our algorithm for 
rotation should reproduce. 

In Figure [H plots of vs. r (top panel) and fl 
vs. r panel (bottom) show that BETHE-hydro per- 
forms well with this test. Simultaneously testing the 
rotational remapping (i 39.7|) and hydrodynamic algo- 
rithms, we employ three remapping regions. The in- 
ner 0.1 cm is Eulerian, zones exterior to 0.2 cm fol- 
low Lagrangian dynamics, and the region in between 
provides a smooth transition between the two do- 
mains. Both and profiles are presented at t — 



0.0, 0.1, 0.2, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, and 1.0 s. As 
designed the angular momentum is conserved to machine 
accuracy. In the top panel (w^ vs. r), it is hard to 
discern any deviation of the simulation (crosses) from 
the analytic solution (solid line). In fact, the maximum 
error as measured by (w^^ana - W0)/niax(v0_ana) ranges 
from ^ 4 X 10^^ near the beginning of the simulation to 
^ 8 X lO"** at the end. O vs. r (bottom panel) shows 
similar accuracy, except for the region near the axis. At 
the center (fiana — ^) reaches a maximum value of 0.2 
rads s~^ at t = 0.55 s. On the whole. Fig. 2] demon- 
strates that this code reproduces the analytic result with 
and without remapping. 

5. GRAVITY 

There are two general strategies for solving Poisson's 
equation for gravity on arbitrary grids. The first is to 
define another regular grid, interpolate the density from 
the hydro grid to this new gravity grid, and use many of 
the standard Poisson solvers for regular grids. However, 
in our experience this approach can lead to numerical 
errors, which may lead to unsatisfactory results. The 
second approach is to solve eq. ([6|) explicitly on the arbi- 
trary hydro grid, eliminating an interpolation step. We 
prefer the latter approach, giving potentials and acceler- 
ations that are more consistent with the hydrodynamics 
on an unstructured mesh. 

In determining the discrete analogs of eq. (O, we 
use the support-operator method. After a bit of alge- 
bra these discrete equations may be expressed as a set 
of linear equations whose unknowns are the discretized 
gravitational potentials. Solving for these potentials is 
akin to solving a matrix equation of the form, Aa; — b, 
and since the solutions change little from timestep to 
timestep, we use an iterative matrix inversion approach 
in which the initial guess is the solution from the previous 
timestep. 

In solving Poisson's equation for gravity, we seek an 
algorithm satisfying the following conditions: 1) the po- 
tential should be defined at zone centers and the gravita- 
tional acceleration should be defined at the nodes; 2) we 
should conserve momentum and energy as accurately as 
possible; and 3) the solver must b e fast. To th is end, we 
emplo y the method described in IShashkov fc S teinberg 
19951) ■ IShashkov fc Steinbergl (jl996f ). and [Morel et al 



19981) for elliptic solvers. Specifically, we fol low the par- 



ticular implementation of iMorel et al.l (|1998f ). 

Applying the support-operator method, we first define 
the divergence operator and, in particular, define (V • 
g)z- Poisson's equation for gravity may be written as an 
equation for the divergence of the acceleration. 



V-g = 4^Gp, 



(48) 



where g, the gravitational acceleration, is the negative of 
the gradient of the potential. 

g^-V*. (49) 
Integrating eq. (|48|) over a finite volume gives 



V • gdV 



g-dS 



AirGpdV , (50) 



IdV JV 

and if the volume of integration, V", is chosen to be that 
of a cell, the RHS of eq. (|50|) becomes —AttGuIz- In 
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discretizing the middle term in eq. ((501 



g-dS 



(51) 



an expression for the divergence of g foUows naturaUy: 



^ e£S(z) 



(52) 



where is the area of edge e, and g^ is the component 
of the acceleration parallel to the direction of the edge's 
area vector. Therefore, substituting eq. (I52p into eq. 
([5D|) suggests a discrete expression for eq. pS)) : 



(53) 



Having defined the discrete divergence of a vector, we 
then use the support-operator method to write the dis- 
crete equivalent of the gradient. Once again, the integral 
identity relating the divergence and gradient operators 
is given by eq. (fT3|) . Applying the definition for the dis- 
crete divergence, the first term on the LHS of eq. (fT3| is 
approximated by 



Jv 



(54) 



where he is the component of H perpendicular to edge 
e, and the RHS of eq. (|13|) is approximated by 



dV 



d5 ~ ^ 



hp 



(55) 



It is already apparent that the discrete potential is de- 
fined for both cell centers ($z) and edges ($e)- This fact 
seems cumbersome. However, we will demonstrate below 
that in the final equations the number of unknowns may 
be reduced by eliminating the cell-centered potentials, 
$2, in favor of the edge-centered potentials, <i>e. Com- 
pleting the discretization of eq. p^ . the second term on 
the LHS is 



(56) 



where Wp is a volumetric weight and Hp and gp are de- 
fined at the nodes. The volumetric weight associated 
with each corner is defined as one quarter of the area 
of the parallelogram created by the sides that define the 
corner. Additionally, the weights are normalized so that 



1. The remaining task is to define the 



dot product Hp ■ gp at each corner when the vector is 
expressed in components of a nonorthogonal basis set: 



H„ — h„ 



(57) 



where Be-i and i-e are basis vectors located at the center 
of the edges on either side of the corner associated with 
node p with orientations perpendicular to the edges, and 
he-i and he are the corresponding magnitudes in this 
basis set. With this basis set, the inner product is 

(58) 



Hp- gp - {Sl,._^Hp, gp) , 



where 



CP 



cost 



cosf 



(59) 



and 9e e-i the angle formed by edges e and e — 1. Using 
eqs. (I57l) - (|59t in eq. (|56)). the second term on the LHS 
of eq. ^ is 

- E -^P ■ dpWpV^ =-^he {{ap + ap+i)ge 

p e 

+Pp9e-1+ Pp9e+l) (60) 

where ap — WpVz/ svc? Op and (3p = 
W^VzCOsOpI SIT? Op. Combining eqs. JM]), ([SS]) 



and (|56p . we have the discrete form of the identity eq. 



[(*e - ^z)heAe - he ((ftp "f Qp+O^e 



E 



+/3p5e-l+/3pge+l)] =0. (61) 

To find the equation for ge on each edge we set the 
corresponding he = 1 and set all others to zero. This 
gives an expression for ge in terms of $e, and the 
other edge-centered accelerations of cell z: 



($e - ^z)Ae - [{ap + ap+i)ge 



[{Up + ap+i)ge + /3p.ge~-l + PpQe+l 



Pp9: 



e+lj 



(62) 



^{^e~^z)Ae. (63) 

Together, the set of equations for the edges of cell z forms 
a matrix equation of the form** 



Q5 = &, 



(64) 



where 6 is a vector based upon the edge values, he — 
Ae{^z — ^e)- Inverting this equation gives an expression 
for the edge accelerations as a function of the cell- and 
edge-centered potentials: 



9e 



'Ae'i^z - $e 



(65) 



After finding the edge acceleration for each zone, these 
values are inserted into eq. ([55)1 . The stencil of po- 
tentials for this equation for each cell includes the cell 
center and all edges. While eq. ([65]) gives the edge ac- 
celerations associated with cell z, there is no guarantee 
that the equivalent equation for a neighboring cell will 
give the same acceleration for the same edge. To ensure 
continuity, we must include an equation which equates 
an edge's acceleration vectors from the neighboring cells: 



Aege,L - AegeM = , 



(66) 



where ge,L and ge.R are the edge accelerations as deter- 
mined by eq. ((^ for the left (L) and right (R) cells. 
While we could have chosen a simpler expression for eq. 
(|66|) . the choice of sign and coefficients ensures that the 
eventual system of linear equations involves a symmetric 
positive-definite matrix. This enables the use of fast and 
standard iterative matrix inversion methods such as the 

On a practical note, all terms involving on the z-axis are 
zero because the areas on the axis of symmetry are zero. Therefore, 
for the zones along the axis, one less equation and variable appear 
in the set of equations for g^. 
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conjugate-gradient method. The stencil for this equation 
involves the edge in question, the cell centers on either 
side, and the rest of the edges associated with these cells. 
At first glance, it seems that we need to invert the system 
of equations that includes the cell-centered equations, eq. 
([65|) . and the edge-centered equations, eq. (|66l) . Instead, 
we use eq. ([55)) to express the cell-centered potentials 
in terms of edge-centered potentials and the RHS of the 
equation: 



(67) 



where Ipe = X^e' ^e^e'e, and = J2e^e' Settee' Ae' = 

J2e' A-e'lpe'- Substituting eq. (pT)) into the expression for 
the edge acceleration, eq. (|65l) . 



9k 



rjk 



AnGm^ + ^ $e^e ( ^V'e - ak. 



Cz 



(68) 



where f?fc = X^e '^fce^e- Finally, substituting eq. (j68p into 
eq. (j66p gives the matrix equation for gravity with only 
edge-centered potential unknowns: 



„R 







Cl 



= -A^GA^Cy-mR+'^mL). (69) 
Qr Cl 

The stencil for this equation involves the edge and all 
other edges associated with the zones on either side of the 
edge in question. Taken together, the equations result in 
a set of linear equations for the edge-centered potentials. 

Since this class of matrix equations is ubiquitous, there 
exists many accurate and fast linear system solvers we 
can exploit. To do so, we recast our discrete form of 
Poisson's equation as 



s, 



(70) 



where $ is a vector of the edge-centered potentials, and 
s is the corresponding vector of source terms. We solve 
the matrix equation using the conjugate-gradient method 
with a multigrid preconditioner. In parti cular, we use the 
algeb raic multigrid package, AMG1R6^ ()Ruge &: StubenI 
flQSl . 

5.1. Gravitational Acceleration 

Unfortunately, this discretization scheme does not ade- 
quately define the accelerations at the nodes where they 
are needed in eq. ([T7|) . Therefore, we use the least- 
squares minimization method to determine the gradi- 
ent on the unstructured mesh. Assuming a linear func- 
tion for the potential in the neighborhood of node i, 
$i(x,?/), we seek to minimize the difference between $fc 
and ^i{xk,yk), where $fc is the neighbor's value and 
^i{xk,yk) is the evaluation of $i(x, y) at the position 
of the neighbor, (xk,yk)- More explicitly, we minimize 
the following equation: 

(71) 



E 



^ki^ik ' 



where the "neighbors" are the nearest cell-centers and 
edges to node i and are denoted by fc e A/'(i), 

= ($, - $fe + $:,,»Aa;,fc -t- ^y^^l:\y^kf , (72) 

^Xik ^ Xk - Xi, Ayik = yk - Vi, and ujf,^ = l/{Axff. + 

^vD- .... _ 

Usually, minimization of eq. (|7ip leads to a matrix 
equation for two unknowns: ^x,i and the gradients 
for <i> in the x and y directions, respectively (^ ^9.2^ . How- 
ever, $i, the potential at the node is not defined and is a 
third unknown with which eq. (|7ip must be minimized. 
Performing the least-squares minimization process with 
respect to the three unknowns, <i>i, and leads 
to the following set of linear equations for each node: 



a$, + 6$^ , + c<^y, = k 



I 



j^yA = rn 



(73) 



where 



k "^ik 



J2k "^tk^^ik 
'Ek ^^k^y^k 



(74) 



^ See i|10.9l for details on performance. 



d = Hk^lk^x^k 
e = 'Ek'^ik^^ik 

9 = Y.k^'ik^^ik^ytk 

h = J2k^ik^yik 

* = Y.k^'ik^^ik^yik 

j = Ek^fk^yfk 

^ = T.k^'ik'^kAx^k 

^ = Hk^'ik'^kAyik 

Inversion of this linear system gives the potential and 
the gradient at the node. Specifically, we use the adjoint 
matrix inversion method to find the inverse matrix and 
the three unknowns, including ^x,i and ^y^i- 

5.2. Gravitational Boundary Conditions 

The outer boundary condition we employ for the Pois- 
son solver is a multipole expansion for the gravitational 
potential at a spherical outer boundary. Assuming no 
material exists outside the calculational domain and that 
the potential asymptotes to zero, the potential at the 
outer boundary, r, is given by 

oo 

0(r) = -G^^ {rrPn{cos9')p{r')dV' , (75) 

n=0 ^ •' 

where Pn{cos9') is the usual Legendre polynomial. In 
discretized form, eq. (j75p becomes 

0(r) = -G^-;^^(r,)"P„(cos0,)m,, (76) 

n—O z 

where Nn is the order at which the multipole expansion 
is truncated. 

5.3. Gravity: ID Spherical Symmetry 

For ID spherically symmetric simulations, calculation 
of gravity can be straightforward using 

G/TT-r) Qncloscd 

9p = ^ , (77) 
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where mp^cnciosod is all of the mass enclosed by point 
p. Instead, we solve Poisson's equation for gravity using 
similar methods to those employed for the 2D gravity. 
By doing this, we remain consistent when comparing our 
ID and 2D results. Furthermore, we show that deriving a 
Poisson solver based upon the support-operator method 
is equivalent to the more traditional form, eq. ([77)1 . 

The inherent symmetries of ID spherical simulations 
reduce eq. to the following form: 



(78) 



Unlike the full 2D version, one may eliminate the edge 
potentials in favor for the cell-centered potentials and the 
accelerations at the edges. Since each edge e will have an 
equation for its value associated with cells z and z + 1, 
this leads to 

$,Ae + V;ge,z = <i>z+lAe + V;+^ ge,z+l • (79) 

Using the continuity of accelerations at the edges, 

9e,z = -9e,z+i , (80) 

the expression for gravitational acceleration at edge e in 
terms of the cell-centered potentials is 



9e,z = rri'^z+i - <i>z)Ae 

Vp 

Substituting this expression into 

Aege^z = -4:TTGm^ , 



(81) 



(82) 



leads to a tridiagonal matrix equation for the cell- 
centered potential, which can be inverted in 0(N) op- 
erations to give the cell-centered potentials. Substitut- 
ing these cell-centered potentials back into cq. ((8T|) gives 
the gravitational acceleration needed for the momentum 
equation (eqs. l^and fTT)) . 

It should be noted that eq. ([82]) can be rewritten as a 
recursion relation for the gravitational acceleration, gp, 
at node p: 



Apgp = Ap^igp^i - 47rGm 



2, interior 



(83) 



where m^. interior is the mass of the cell interior to point 
p. With the boundary condition gi = and then recur- 
sively solving eq. ([55)1 for gp from the center outward, 
the expression for the gravitational acceleration simpli- 
fies to eq. (177)) . This is the traditional form by which 
the gravitational acceleration is calculated for ID spher- 
ically symmetric simulations and for which the potential 
is usually not referenced. The beauty of our approach is 
that ID and 2D simulations are consistent, and that it 
self-consistently gives the gravitational potential. 

5.4. Conservation of Energy with Gravity 

For self-gravitating hydrodynamic systems, the total 
energy is 

E = £ + K + U , (84) 

where £ is the total internal energy, K is the total kinetic 
energy, and U is the total potential energy: 



U : 



1 



p<^dV, 



(85) 



An alternative form is possible with the substitution of 
eq. pSj) into eq. ([55)1 and integration by parts: 

1 



U ■ 



8nG 

The integral form of eq 
E^ I {pe 



g^-dS+ / \g\'dV 
s Jv 



IV 2" 2 
and total energy conservation is given by 



dE d 



1 



dt dt.ly^P'^r''^ 



1 



-p^)dV = 0. 



(86) 



(87) 



(88) 



Since the solution to Poisson's equation for gravity, using 
Green's function, is 



$(f) = -G 



p{x')dV' 



(89) 



the time derivative of the total potential energy may take 
three forms: 



V dt 



dt 2 



d{p<P) 
dt 



dV . 



(90) 



Note that this is true only for the total potential energy 



and does not imply = = Using eq. ^ 

and the Eulerian form of the hydrodynamics equations 
(see appendix [A]) . eq. (l88|) becomes a surface integral: 



P + \pv^ 



p<^)v ■ dS , 



(91) 



which simply states that energy is conserved in the ab- 
sence of a flux of energy through the bounding surface. 

The discrete analogs of £ and K are trivially obtained 
using £ — ^^rnzSz and K = ^X^p^p'^p; respectively. 
The discrete total potential energy may take two forms, 
the analogs of eqs. ((85|) and ((86)) : 



U = 



2 



1 

SttG 



rriz^z 

J2pb ^Pb9pb ' 



Spb + J2p 9p • 9pVp 



(92) 



where pb indicates outer-boundary values. Upon trying 
both forms, we get similar results. Therefore, we use 
the simpler form involving the cell-center potential: U = 

\Y.z^^^^- 

Because our discrete hydrodynamics equations includ- 
ing gravity do not give strict energy conservation, we use 
the core-collapse simulation f i|10.9p to gauge how well 
energy is conserved. Rather than measuring the relative 
error in total energy by AE/Ei-d, where AE = En — E^-d, 
E-cci is the initial total energy, and £"„ is the total energy 
at timestep n, we use Ai?/|[/„|. For stellar profiles, the 
kinetic energy is small and the internal and gravitational 
energies are nearly equal in magnitude, but opposite in 
sign. Since the total energy is roughly zero in compar- 
ison to the primary constituents, £ and U , of the total 
energy, we use the gravitational potential energy as a 
reference. For example, the internal and gravitational 
energies start at ~ 4 x 10^^ erg and reach 1 x 10^'^ 
ergs at the end of the run. However, the total energy is a 
small fraction of these energies initially, 5 x 10^°, and 
after core bounce, ~ 1.5 x 10^^. Hence, we measure the 
relative error in total energy as Ai?/|[/„|. 
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The total energy for the ID and 2D core-coUapse sim- 
ulations evolve similarly and is conserved quite well. For 
all times except a few milliseconds around bounce, total 
energy is conserved better than ~1 xlO^^. During col- 
lapse, from to 147 ms, the total energy deviates by only 
~3 xlO"''. Over a span of 5 ms around bounce, t = 148 
ms, the total energy changes by ~2.2 xlO~^. For 100s 
of milliseconds afterward, the relative error in the total 
energy reduces to ~7 x 10""^. Hence, for all but about 5 
ms of a simulation that lasts many hundreds of millisec- 
onds the relative error in total energy is better than ~1 
xlO-3. 

5.5. Tests of Gravity 

We include here several assessments of the ID and 2D 
Poisson solvers. First, we calculate, using the butterfly 
mesh, the potential for a homogeneous sphere and com- 
pare to the analytic solution. This tests overall accuracy 
and the ability of our solver to give spherically symmetric 
potentials when a non-spherical mesh is employed. Sim- 
ilarly, we substantiate the algorithm's ability to produce 
aspherical potentials of homogeneous oblate spheroids. 
Then, we verify that the hydrodynamics and gravity 
solvers give accurate results for stars in hydrostatic equi- 
librium, and for a dynami cal problem, the Goldrei ch- 
Weber self-similar collapse (jGoldreich fc Web"erl|198 tf). 

Figure [5] compares the simulated and analytic solutions 
for a homogeneous sphere, having density po = 1 g cm"'^, 
and maximum radius = 1 cm. The analytic potential 
inside the sphere is <i>ana — Gpo^T^ir^ H- — 3r^), and 
the top panel of Fig. [5] shows the relative difference of 
the analytic potential and the numerical solution, <I>z, as 
a function of radius. Results for four resolutions of the 
butterfly mesh are shown: 2550 cells with an effective 
radial resolution Ar 0.02 cm (blue); 8750 cells with 
Ar ~ 0.01 cm (green); 15,200 cells with Ar ~ 0.006 cm 
(yellow); and 35,000 cells with Ar ~ 0.005 cm (red). Two 
facts are obvious: 1) the solutions are accurate at a level 
of - 3 X lO"'^ for Ar/ra ^ 0.005 and - 3 x 10"'' for 
Ar/va ^ 0.02; and 2) the degree to which the solution 
is spherically symmetric is similarly a few times 10^^ 
for Ar/ra ^ 0.005 and a few times 10"'' for Ar/ra ~ 
0.02. Plotting the minimum error as a function of the 
effective radial resolution, the bottom panel of Fig. [5] 
verifies that the 2D Poisson solver convergences with 2"^^- 
order accuracy (solid line). 

Next, we calculate the aspherical potential for an 
oblate spheroid . The homogeneous oblate spheroid has 
an elliptic meridional cross section, and the minor- and 
major-axes of the ellipse are rf, and r^, where is the 
equatorial radius. Thus, the eccentricity of the spheroid 
is e = -^/l — (rb/ra)'^- Given a uniform density, po, the 
potential for a spheroid is 

$ana(r, z) = -nGpo {IbttI - [flir^ + agz^]) , (93) 

where, for oblate spheroids. 



fli = 



as = 2 



sin ^ e 



-(1-e^) 



2n1/2 



(l-e2)V2 



(l-e2)-i/2 



sin ^ e"' 



(l-e2)i/2 



(94) 



(95) 



and Ibt — 2ai + 03(1 — e^) l)ChandrasekhaH[l969f) . 



Numerical results with e = 0.8 are shown in Fig. [6l As 
with the sphere, po = 1 g cm^"^ and ra — 1 cm. However, 
for the given eccentricity, the polar-axis radius, r^, is 0.6 
cm. Once again, the grid is a butterfly mesh, but this 
time the outer boundary follows the ellipse defining the 
surface of the spheroid. The top panel of Fig. [6] presents 
the spheroid's potential, and the degree of accuracy is 
presented in the bottom plot. With iVccii = 35, 000 and 
Ar/ra ~ 0.005, the relative error in the potential ranges 
from ^ 2 X 10~® near the outer boundary to ~ 3 x 10~^ 
in interior regions. Conspicuous are features in the rel- 
ative error that track abrupt grid orientation changes in 
the mesh. Fortunately, these features have magnitudes 
smaller than or similar to the relative error in the local 
region. The relative error in the gravitational acceler- 
ation magnitude, {\gp\ - Iganal/maxd^anal) ranges from 
~— 10"* to ^10"'*, where ^ana is the analytic accelera- 
tion and max(|gana|) is the maximum magnitude on the 
grid. Typical errors in the acceleration direction range 
from ~— 10"''^ to ~10"* radians with rare deviations as 
large as ~10~'^ radians near the axis and abrupt grid 
orientation changes. 

In Fig. [71 we demonstrate that the ALE algorithm 
in combination with our gravity solver produces reason- 
ably accurate hydrostatic equilibria. The grid is the but- 
terfly mesh with 8750 zones, and the initial model is a 
Lane-Emden polytrope with j — 5/3, M — 1 Mq, and 
R = 2.9 X 10^° cm. Crosses show the density profile at 
t = 1 X 10* s, while the solid line shows the maximum 
density as a function of time and that the star pulsates. 
These oscillations result from the slight difference be- 
tween an analytic hydrostatic equilibrium structure and 
a discretized hydrostatic equilibrium structure, and with 
increasing resolution, they decrease in magnitude. Inter- 
estingly, the oscillations continue for many cycles with 
very little attenuation. 

Simulating the Goldrei ch- Weber self-similar collapse 
(jGoldreich fc WebedllMih in ID and 2D is a good test of 
dynamic simulations including gravity. The analytic pro- 
file is similar to a Lane-Emden polytrope with 7 — 4/3, 
and in fact, we use the gamma-law EOS with 7 — 4/3. 
While the Lane-Emden polytropes are assumed to be in 
hydrostatic equilibrium, the Goldreich- Weber self-similar 
collapse has a homologous velocity profile and a self- 
similar density profile. The physical dimensions have 
been scaled so that M = 1.3 Mq, the initial central den- 
sity is 10^" g cm~"^, and the maximum radius of the pro- 
file is 1.66 X 10* cm. For the ID simulation, we initiate 
the grid with 200 evenly-spaced zones, and for 2D, a but- 
terfiy mesh with 35,000 zones (effectively with 200 radial 
zones) is used to initiate the grid. Subsequent evolution 
for both simulations uses the Lagrangian configuration. 

Figure [8] shows snapshots of density vs. radius for ID 
(top panel) and 2D (bottom panel) simulations at t —0, 
20, 40, 60, 80, 100, 120, and 130 ms. Both plots indicate 
that the simulations (crosses) track the analytic solution 
(solid lines). Quantitatively, we measure a relative dif- 
ference, {pz - Panay(max(pana)), for the reported times, 
where pz are the simulated cell-center densities and pana 
are the analytic values. Consistently, the largest devi- 
ations are at the center and Table [2] gives these values 
for the ID simulation (2"'' column) and the 2D simula- 
tion (3'''' column). The relative differences range from 
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—9.2 X 10 ^ at ms to -3.6 X 10^2 at 130 ms for the ID viscosity tensor is a combination of a scalar coefficient, 
simulation and —1.2 x 10^^ at ms to —6.3 x 10^^ at 



130 ms for the 2D simulation. At all times, the departure 
from spherical symmetry is no more than 1 x 10^*. 

6. HYDRODYNAMIC BOUNDARY CONDITIONS 

Boundary conditions are implemented in one of two 
ways. Either an external pressure is specified or the ve- 
locities at the nodes are fixed. For external pressures, 
ghost cells are defined that have no true volume or mass 
associated with them. Their only function is to apply an 
external force to the boundary cells equal to an exter- 
nal pressure times the boundary surface area. Specifying 
nodal velocities on the boundary accomplishes the same 
task. In this case, there is an implied external force and 
pressure. In practice, generic dynamical boundaries use 
the external pressure boundary condition. On the other 
hand, pistons, reflecting walls, and the azimuthal axis 
have the velocity perpendicular to the boundary speci- 
fied, while the parallel component executes unhindered 
hydrodynamic motions. 

7. ARTIFICIAL VISCOSITY 

To resolve shocks over just a few zones, we include arti- 
ficial viscosity in the equations of hydrodynamics. For ID 
simulations, we add a v iscous-like term to the pressure 
(jVon Neuman n fc Richt mver 1950). We denote this vis- 
cous pressure by q. The original realization ofq employed 
one term proportional to (Au)^, where Av is the differ- 
ence in velocity from one zone to the next. While this 
form adequately resolved shocks, unphysical oscillations 
were observed in the post-shock flow. A second term, 
linear in Au, wa s then added th at effectively damped 
these oscillations (lLandshofilll955[ ) . Another form which 
we employ for our ID simulations is 



plC2^^\Av\ 




(Aw)2 + c2c2 



lAt'l, (96) 




HOC p \ C2 



(97) 

and the gradient of the velocity tensor, Vu. Therefore, 
references to the parameters ci and C2 in ID and 2D 
refer to the same linear or quadratic dependence on \Av\. 
With this assumption, the momentum equation, eq. ([2]), 
becomes 



dv ^ ^ , ^ 

p— = -pv$ - VP + V • (^iVv) 

at 

and the energy equation, eq. ([3]), is 

de ^ ^ - 

p— = -PV • v + p.(Vv) : (Vw) 
dt 



(98) 



(99) 



where (Vu) : (Vw) = (Vt7).y (Vi/)^ , using the normal 
Einstein summatio n convention. With t h e me thod of 
support operators, I Campbell fc Shashkovl (|2001h derive 
discrete forms of these equations. Analogous to the gra- 
dient of the pressure term, the discrete artificial viscosity 
term in the momentum equation becomes a summation 
of corner forces. 



V • (/iVt?) 



E 



FP 



(100) 



and the corresponding term in the energy equation is the 
usual force-dot-velocity summation: 



/l(V^/) : (Vv) 



E 



/p.visc ■ 



(101) 



Therefore, implementation of the artificial viscosity 
scheme is straightforward and similar to that for the pres- 
sure forces. 

In practice, this artificial viscosity scheme is formu- 
lated for Cartesian coordinates, and to obtain the equiv- 
alent force for cylindrical coordinates, we multiply the 
Cartesian subcell force by 27rrp. In general, this works 
quite well, but sacrifices strict momentum conservation 
(See SU). 



where Cg is the sound speed, ci is the parameter associ- 
ated with the linear term, and C2 is the parameter asso- 
ciated with the quadratic term. This form has the ap- 
pealing attribute that it is motivated by the expression 
for the shock- jump condition for pressure in an ideal gas 
fWilkinslll980IV 

For 2D simulations, the artificial viscosity scheme we 
have settle d upon is the tensor a rtifici al viscosity al- 
gorithm of iCampbell fc Shashkovl (l200lh . Here, we do 
not re-derive the artificial viscosity scheme, but high- 
light some of its salient features and practical imple- 
mentations. The useful feature of the tensor algorithm 
is its ability to calculate artificial viscosity on an arbi- 
trary grid, while suppressing artificial grid buckling when 
the flow is not aligned with the grid. In the strictest 
sense, the tensor artificial viscosity does not employ 
the simple viscous pressure described above. Instead, 
ICampbell fc Shashkovl (|200lh assume that the artificial 



8. SUBCELL PRESSURES: ELIMINATING HOURGLASS 
GRID DISTORTION 

A problem that can plague Lagrangian codes using 
cells w ith 4 or more sides is the u nphysical hourglass 
mode (jCaramana fc Shashkovl [l998l ). To suppress this 
problem, we employ a modified version of the sub - 
cell pressure algorithm of lCaramana fc Shashkovl (|1998f) . 
For all physical modes (translation, shear, and exten- 
sion/contraction) the divergences of the velocity on the 
subcell and cell levels are equal. 



(V • v)^ = (V • v). 



(102) 



Hence, as long as the calculation is initiated such that the 
subcell densities of a cell are equal to the cell-averaged 
density, then they should remain equal at all subsequent 
times. Any deviation of the subcell density from the cell 
density, 6pp — Pp — pz, is a direct result of the hourglass 
mode. The scheme of lCaramana fc Shashkovl (119981 ) uses 
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this deviation in subcell density to define a subcell pres- 
sure that is related to the deviation in subcell and cell 
density, by 

P^^P,+clSp; (103) 
with the corresponding deviation in pressure: SP^ — 

Application of these pressures as subcell forces coun- 
teracts the hourglass distortion. However, the pressure 
throughout the cell is no longer uniform. As a result, the 
subcell pressures exert forces on the cell centers and mid- 
edge points. These locations in the grid are not subject 
to physical forces as are the nodes. Instead, their move- 
ment is tied to the motions of the nodes. To conserve 
momentum, the forces acting on these enslaved points 
must be redistributed to the dynamical nodes. Though 
the choice of redistributi on is not uniqu e, we follow the 
procedure established by ICaramana fc Shashkovl p998f ) 
for the form of the subcell force (for the definitions of Sp 
and flp, see S|3]and Fig. [S]): 

f; = sp^s; + \ [{5p; - 5p;^,)d; + {5p;_, - 5p;)a;_,] . 

This particular scheme works quite well for Lagrangian 
calculations. However, it can be incompatible with sub- 
cell remapping algorithms. The hourglass suppression 
scheme described above assumes that only hourglass mo- 
tions result in nonzero values for 5pp. Even in calcula- 
tions that are completely free of hourglass motions, the 
subcell remapping scheme can and will produce subcell 
densities within a cell that are different from the cell den- 
sity. Since this difference of subcell and cell densities has 
nothing to do with hourglass motions, any subcell forces 
that arise introduce spurious motions that don't correct 
the hourglass motions. 

Investigating severa l alternative schemes, we settle d on 
a modification of the ICaramana fc Shashkovl (|1998D ap- 
proach. After each remap, we define a tracer density 
whose only purpose is to track the relative changes in the 
subcell and cell volumes during the subsequent hydrody- 
namic solve. For convenience, and to keep the magnitude 
of the subcell pressures about right, we set the subcell 
tracer density equal to the cell density. After replacing 
5pp with the difference between the subcell tracer density 
and the cell d ensity, 6p^j.j., we procee d with the scheme 
prescribed by ICaramana fc Shashkovl (Il998f). We have 
found that the subcell forces of ICaramana fc Shashkovl 
()1998[ ) significantly resist hourglass motions only after 
6pp achieves significant magnitude. For our new hour- 
glass scheme, if a remap is implemented after each hydro- 
dynamic solve, Spp ^^ does not have a chance to achieve 
large values, and hence, the subsequent subcell forces 
are not very resistant to hourglass distortions. However, 
for many simulations, we have found that it is not nec- 
essary to remap after every Lagrangian hydrodynamic 
solve. Rather, the remap may be performed after N 
timesteps. In these circumstances, the subcell tracer den- 
sity is allowed to evolve continuously as determined by 
the hydrodynamic equations. With large N, Spp does 
develop significant amplitude and the subcell forces be- 
come effective in suppressing the hourglass modes. In 
fact, in the limit that ^ 1, this scheme becomes the 



hourg lass suppression scheme of ICaramana fc Shashkovl 
(fT99l . 

In practice, we multiply this subcell force by a scal- 
ing factor. To determine the appropriate magnitude of 
this scaling factor, we executed many test problems and 
found that problems that involved the perturbation of 
hydrostatic equilibrium provide a good test of the robust- 
ness of the hourglass fix. During the testing protocol, we 
ran a simulation for at least 200,000 timesteps, ensuring 
that the hourglass fix remains robust for long calcula- 
tions. For Lagrangian calculations, a scaling factor ~ 1 
produced reasonable results. For runs in which remap- 
ping occurred after N timesteps, we found that larger 
values of the scaling factor were required. However, for 
very large values of N, such as 64 or greater, such large 
values compromised small structures in the flow. Hence, 
we settled on scaling factors near 1 for large N. For small 
N, we found that the difference in tracer density and the 
cell-centered density had little time to build to signif- 
icant amplitudes, so larger values of the scaling factor 
are required, but anything above 4 produced noticeable 
problems in flows with many timesteps. In summary, the 
scaling factor should be adjusted to be ~ 1 for large N 
and ~ 4 for small N. 

8.1. Tests of the Hourglass Elimination Algorithm 

To demonstrate the need for an hourglass suppression 
scheme, and that our algorithm to address the hour- 
glass distortions works, we show in Fig. [9] results for the 
Goldreich- Weber self-similar collapse in 2D fi iS.Sp . with 
and without the hourglass suppression (Sj8]). These re- 
sults represent a Lagrangian simulation of the self-similar 
collapse at t — 118 ms. On the left ("r < 0") the grid 
clearly shows problematic grid buckling when the hour- 
glass suppression is turned off. The grid on the right 
("r > 0") shows the elimination of grid buckling with 
the use of the hourglass suppression scheme and a scal- 
ing factor of 2.0. 

FigurefTUl illustrates the problem of using the hourg lass 
suppression scheme of ICaramana fc S hashkov (1998) in 
combination with the subcell remapping algorithm (fj9]). 
All three panels show the results of the single-mode 
Rayleigh- Taylor instability at t = 12.75 s. For the left 
panel, no hourglass suppression is employed. This panel 
represents our control. While the hourglass instability 
does not cause serious problems for the calculation, one 
can see evidence of slight hourglass patterns at the scale 
of the grid resolution. In particular, two blobs near the 
top and on the edges form distinct patterns. The cen- 
tral panel demonstrates the problem of using the subcell 
remapping scheme and the s ubcell pressure method of 
ICaramana fc Shashkovl ()1998f ). On the other hand, the 
right panel demonstrates that the modifled subpressure 
scheme we have developed suppresses the hourglass dis- 
tortions, while preserving the expected flow. 

9. REMAPPING 

Fluid flows with large vorticity quickly tangle a La- 
grangian mesh, presenting severe problems for La- 
grangian hydrodynamic codes. To avoid entangled grids, 
it is common to remap the state variables after each La- 
grangian hydrodynamic advance to another less tangled 
mesh. Even flows that exhibit very little vorticity, but 
extreme compression or expansion along a particular di- 
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rection, leading to very skewed cells, can limit the ac- 
curacy of Lagrangian calculations. In this case as well, 
one can employ a remapping scheme which is designed 
to minimize the calculational error. 

To remap, a new grid must be established. One choice 
is to remap to the original grid, thereby effectively solv- 
ing the hydrodynamics equations on an Eulerian mesh. 
Another is to use a reference Jacobian matrix rezone 
strategy, which establishes a new mesh close to the La- 
grangian mesh, while reducing numerical error by reduc- 
ing unnecessarily skewed zones. The arbitrary nature 
of the remapping algorithm even allows the freedom of 
choosing one rezoning scheme in one sector of the grid 
and another rezoning scheme in other sectors. For exam- 
ple, one could remap in an Eulerian fashion in one sector, 
let the grid move in a Lagrangian manner in another, and 
in the intervening region smoothly match these two re- 
gions. 

Generally, there are two options for remapping 
schemes. One can remap from one arbitrary grid to an- 
other completely unrelated arbitrary grid. For these un- 
related grids, one needs to determine the overlap of the 
zones of the first grid with the zones of the second. In 
general, this can be a very cumbersome and expensive 
process. As a result most ALE codes, remap to an arbi- 
trary grid that is not too different from the first. Specif- 
ically, the usual stipulation employed is that the face of 
each cell does not traverse more than one cell during a 
timestep, and that the connectivity among nodes, faces, 
edges, and cells remains the same from the first grid to 
the remapped grid. The regions swept by the faces then 
contain the mass, momentum, or energy which is added 
to the cell on one side of the face and subtracted from 
the cell on the other side. 

It is in this context t hat we use the swcpt-region 
remapping algorithms of ILoubere fc Shas hkov (200. 



In order to remap the primary quantities using the 
subcell remapping algorithm, mass, momenta, internal 
energy, and kinetic energy must be defined at the subcell 
level. In the Lagrangian stage, the subcell mass, m^, is 
already defined. The subcell density which is a conse- 
quence of the change in the subcell volume (Vp) during 
the Lagrangian hydro step is 

Pl-y^- (105) 

Unfortunately, there is no equivalent quantity for the 
subcell internal energy. Instead, we define the internal 
energy density for each zone, = SzPzi fit a linear func- 
tion for the internal energy density, — + (Ve)^ • (x — 
Xc), and integrate this function over the volume of each 
subcell to determine the total subcell internal energy and, 
consequently, the average subcell internal energy density: 

E.^ - edV and = ^ • (106) 

Performing the integration of eq. p06p involves a vol- 
ume integral and volume integrals weighted by x and y 
in Cartesian coordinates or r and z in cylindrical coordi- 
nates (see appendix [B1 for formulae calculating these dis- 
crete integrals). Note, that, by construction, these newly 
defined internal energies satisfy conservation of internal 
energy for each zone, m^e^ = J2pesiz) ^p- 

Similarly, there are no readily defined subcell-centered 
averages for the momenta, flp, or velocities, Up. However, 
they should be related to one another by 

fi; - m^ii; . (107) 

If the gathering stage is to preserve momentum conser- 
vation, then the following equality must hold: 

5Z "^p"p = H ^p'"p ■ 
peS{z) pes{z) 



(108) 



Margolin fc Shaskovl (l2003l . l2004l) . iKucharik et al.l By inspection, it might seem natural to set Up 



(|2003f) . and ILoubere etlir (|2006D " This remapping 



scheme may be described by four stages: 

1. The first is the gathering stage. This is a stage 
in which subcell quantities of the mass, momenta, 
and energies are defined in preparation for the bulk 
of the remapping process. 

2. The second is to remap the subcell quantities from 
the Lagrangian grid to the new rezoned grid using 
the swept-region approximation. In doing so, the 
remapping algorithm remains 2nd-order accurate, 
but avoids time intensive routines to calculate the 
overlap regions of the old and new grids. 

3. The third is to repair the subcell densities. Be- 
cause exact spatial integration is avoided with the 
swept-region approximation, local bounds of sub- 
cell densities may be violated. Therefore, a re- 
pair algorithm which redistributes mass, momen- 
tum, and energies to preserve the local bounds is 
implemented. 

4. Finally, there is the scattering stage, in which the 
primary variables of the hydrodynamics algorithm 
are recovered for the new rezoned grid. 

9.L Gathering Stage 



Ho wever, we follow the more a ccurate suggestion made 
by ILoubere &: Shashkovl (|2005( l to define the subcell- 
averaged velocity as follows: 

^ ^ + + W + (109) 

where is obtained by averaging over all nodal velocities 
associated with zone z, and and Vp-1/2 are edge- 



averaged velocities given by t/p+i/2 = 1/2 {vp + Vp+i) and 
Vp-1/2 = 1/2 {vp + Vp_i). Substituting eq. (|109p into eq. 
IjlOSp and rearranging terms, we get an expression for the 
subcell velocity that depends on known subcell masses 
and nodal velocities: 



Vp-l 



'^{Avp, -Vp>+i-Vp>-i). (110) 



p'£S(z) 

We now rewrite eq. (jllOp in a form that obviously makes 
it easy to write the gathering operation in matrix form: 



u- 



E 

p'£S{z) 



8m, 



(111) 
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If IJz is a vector of one component of all velocities asso- 
ciated with cell z, and Uf*^ is the equivalent for subcell 
velocities, then 

Ufc-^M.U,, (112) 

where the matrix, M^, for each zone is given by the co- 
efficients in eq. (jllll) . 

To conserve total energy in the remap stage, we define 
the subcell kinetic energy and follow the same gather 
procedure as for the velocity by simply replacing velocity 
with the specific kinetic energy: 



(113) 



We ensure conservation of total kinetic energy for a cell: 



E 



iripkp — 



E 



(114) 



where is the subcell-averaged specific kinetic energy. 
Then, the transformation from nodal specific kinetic en- 
ergies to subcell-averaged specific kinetic energies is given 

by 

kfc- ^ M,k, , (115) 

where kf '-^ is a vector of the subcell-averaged specific ki- 
netic energies, is a vector of the nodal specific kinetic 
energies, and is the same transformation matrix used 
for the velocity transformations. Having found the spe- 
cific kinetic energies, the calculation of the kinetic energy 



for each subcell is then straightforward: 



Ikt. 



Upon completion of the gathering stage, each relevant 
quantity, m^, flp, E^, and K^, is expressed as a fun- 
damentally conserved quantity for each subcell. From 
these conserved quantities and the subcell volumes, we 
have the corresponding densities, which, as is explained 
in ^9.2\ are important components of the remapping pro- 
cess. 

9.2. Swept-edge Remap 

Having gathered all relevant subcell quantities, we be- 
gin the bulk of the remapping process. Since all vari- 
ables are now expressed in terms of a conserved quantity 
(subcell mass, energy, and momentum) and a density 
(mass density, energy density, and momentum density), 
for clarity of exposition, we focus on representative vari- 
ables, subcell mass and density, to explain the generic 
remapping algorithm. When there are differences in the 
algorithm for the other variables, we note them. 

We have in s titute d the remapping algorithm of 
iKucharik etal\ ()2003f ). To avoid the expensive process 
of finding the overlap of the Lagrangian grid with the 
rezoned grid, this algorithm modifies the mass of each 
subcell with a swept-edge approximation for the rezon- 
ing process. As long as the connectivity and neighbor- 
hood for each cell remain the same throughout the cal- 
culation, time spent in connectivity overhead is greatly 
reduced. The objective of the remapping algorithm is 
then to find the amount of mass, Sirie, that each edge 
effectively sweeps up due to the rezoning process and to 
add or subtract this change in mass to the old subcell 
mass to obtain the new subcell mass, mf;: 



(116) 



Of course, the change in mass is obtained by integration 
of a density function over the volume of the swept region: 



pdV. 



(117) 



The density function, p{x, y), may take on any functional 
form. To maintain second-order accuracy of the algo- 
rithm we use a linear density function; 



p = pf + (Vp)f .(f-fP) 



(118) 



where is the average of the nodal positions defining 
the subcell. 

In determining the gradient, (Vp)f , we employ one of 
two standard methods for calculating the gradient using 
the densities from subcell {z,p} and its neighbors. The 
first method uses Green's theorem to rewrite a bounded 
volume integral as a boundary integral around the same 
volume. We begin with a definition of the average gra- 
dient of p for a given region. For each component, the 
average gradient is 



dx / ~ V JV f^a:"!' — V 

- vIvPydV - - 



TdV 



■dp 
, 9v 



pdy 



V §dv ' 



(119) 



where p^ and py are the gradients in the x- and y- 
directions, respectively. The region over which the in- 
tegrals are evaluated is defined by segments connecting 
the centers of the neighboring subcells. Hence, the dis- 
crete form of the average gradient is 



dx 



dpY _ 



(120) 



While this method works reasonably well, an unfortunate 
consequence of the integration and division by volume is 
that the value of the gradient is directly influenced by 
the shape and size of the cell. 

An alternative method for determining gradients, 
which is not as easily influenced by the shape of the 
mesh, is a least-squares procedure. For simplicity of no- 
tation, we describe this method in the context of cells 
rather than subcells. Again, assuming a linear form 
for pz{x,y), we seek to minimize the difference between 
Pc and pz{xc,yc), where pc is the neighbor's value and 
Pz{xc,yc) is the extrapolation of cell z's linear function 
to the position of the neighbor, (xc,yc)- More explicitly, 
we wish to minimize the following equation: 



CZ ZC 1 



(121) 



ceA/'(2) 



where the set of neighbors for cell z is denoted by c € 

AA(z), 

Ej^ = (-Ape. + Px.zAXcz + Py,z^yczf , (122) 

Px.z and py z are the x and y components of the gradient, 
Apc2 = Pc-Pz, Axc^ = Xc-Xz, Ayc z = y c-yz, and uj'^^ = 
l/{Axl^ + Ay^^). Minimizing eq. (|12ip with respect to 
the two unknowns, leads to the following set of linear 
equations: 



eGS(z) 



apx + bpy = d 
bpx + cpy = e , 



(123) 
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for each subcell, where 



a 


— J2ceAr{z) 


^Iz^xl, 


b 




^Iz^^cz^Vcz 


c 


= J2ceAf{z) 


^cz^Vcz 


d 




UJ^^Apcz^Xcz 


e 


— J2cGAr{z) 


<^cz^PczAycz 



(124) 



Solving this hnear system with Cramer's rule then gives 
the least-squares gradients. 

After calculating the gradient we ensure mo notonicity 
using the Barth-Jespersen limiter (|Barthlll997f ). The gra- 
dient is limited by a scalar, that has a range between 
and 1: 

p,{x,y)=p, + <i>,{Wp),-{x~x,). (125) 

First, we determine the minimum and maximum value 
of p among pz and its neighbors pc'. 



Pz 



min(p^,Pc 



pr''-max(p^,Pe). 
We satisfy the requirement that 

min^ ^ / \ ^ max 

Pz <p[x,v)<Pz ■ 
This is accomplished in the following way. 

= min($^) , 



(126) 
(127) 

(128) 

(129) 



where is a limiter associated with each node of the 
cell and is given by 



min ( 1, 



Pz 



min ( 1, 
1 



for pixn,yn) - Pz > 

for piXn,yn) - Pz <0 

for p{xn,yn) - Pz ^ , 
(130) 

where /5(x„, y„) is the unlimited linear hmction evaluated 
at the nodes. 

With the linear function p defined for each subcell, we 
return to the task of determining the swept mass, Srue- 
The specific linear function p for density used in eq. pi7p 
depends upon which subcell the edge encroaches upon. 
In fact, it is quite possible that the volume swept by the 
edge intersects more than one immediate neighbor of the 
subcell in question. Of course, an accurate swept-edge 
integration scheme would take i nto account all spatia l 
functions in the relevant subcells. iKucharik et al.l (|2003l ). 
and references therein, have noted that such an accurate 
scheme may be cumbersome and computational expen- 
sive. Instead, they suggest an approximate swept region 
integration in which only the subcells to the left and right 
of the oriented edge need matter in determining the spa- 
tial function of density used for integration. To this end, 
an oriented volume integral is defined: 



6Ve 



xdy . 



(131) 



Whether this volume integral is negative or positive, the 
density function is taken from either the left or right 
subcell: 

^_jPr, SVe>0, 
P-\ PI, SVe<0 



(132) 



9.3. Repair 

While the above scheme is second-order accurate and 
mass conserving, the approximations made in the swept- 
edge algorithm can violate local bounds. In practice, the 
values of the remapped densities should be bounded by 
the neighborhood values of the original grid. Computa- 
tionally expensive schemes that find the overlap regions 
among the old and new grid can be made to preserve the 
bounds. However, the local bounds can be violated in 
swept-edge remappin g. Using the linearity -and-bound- 
preservmg method of IKucharik et all ('2003'), we conser- 
vatively repair quantities by redistributing mass, energy, 
momenta, and number to neighbors satisfying the local 
bounds of the quantities on the previous Lagrangian grid. 
However, the order in which one processes the cells and 
subcells influences the specific values that emerge from 
the repair process. Therefore, we employ the order in- 
dependent scheme of iLoubere et al.l ()2006f ) . These repair 
schemes are designed to repair subcell densities. For ve- 
locities and composition we wish to preserve the bounds 
of V and Xi. Therefore, we have appropriately adjusted 
the scheme to conservatively repair momenta and parti- 
cle number by preserving the local bounds of v and Xi, 
respectively. 

9.4. Scattering 

The final remap step is to recover the primary variables 
for the next hydrodynamic step. The new cell-centered 
and node-centered masses are 



E 

pesiz) 



and 



2es(p) 



(133) 



(134) 



respectively. Consequently, the subcell and cell-averaged 
densities are 

p; = ^ (135) 



and 



rhz 

'Vz 



(136) 



respectively. 

Recovering the velocities at the nodes is slightly more 
involved. First, the new subcell-averaged velocity is de- 
fined using the remapped momenta and masses: 



Pi 

mi' 



(137) 



p 



Then, for each cell we invert the matrix equation that 
transforms node velocities into subcell velocities: 



(138) 



In this scheme, each cell provides its own velocity for a 
node, which may be different from another cell's value 
for the same node. Therefore, for each node velocity, we 
average node velocities resulting from the inversion of eq. 
(PSll : 



zes{p) 
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Finally, we recover the specific internal energy. Com- 
paring the sum of the old internal energy and kinetic 
energy with the sum of remapped internal energy and a 
kinetic energy defined by the remapped node velocities 
will not necessarily ensure conservation of total energy. 
Therefore, the old kinetic energy is remapped along with 
the node velocities. The discrepancy between these two 
kinetic energy representations is then added to the inter- 
nal energy, ensuring conservation of total energy during 
the remap stage: 



£r. 



peS(z) 





(140) 



The new cell-centered specific internal energy is then ob- 
tained using the modified internal energy, £z, and the 
remapped mass, m^: 



(141) 



This ensures conservation of total energy at the 
expense of consistent remapping of internal energy. 
Whichever is desirable depends upon the problem. For 
example, in flows with large kinetic energies and small in- 
ternal energies, the discrepancy in kinetic energy remap- 
ping could substantially alter or even dominate the inter- 
nal energy. Consequently, we include a flag in BETHE- 
hydro that determines whether the remapped kinetic en- 
ergy differences are added to the internal energy. 

9.5. Remapping Tests 

Figure [11] demonstrates the basic effectiveness of the 
remapping algorithm. We remap a ID step function in 
density many times. From a; = to x = 1/2 cm, the 
density is p = 2.5 g cm~^, and from a; = 1/2 to x = 1 
cm, the density is p = 1.5 g cm~'^. This test has 50 cells 
and Nnodes = 51 nodes. Indexing each node by i, the 
positions of the nodes are 



Xi = {I - a) 



where 



i - 1 



i- 1 



■^nodes 1 



1 . / n 

a = — sin 47r 



, (142) 



(143) 



and the remapping step is n. As a result, the grid 
completes two full cycles in this remapping test. Fig- 
ure [TT] shows the result of this remapping test where 
nmax = 800. The top panel displays the density pro- 
file for remapping steps to Umax/'^, while the bottom 
panel shows the profile for steps nmax/^ to Umax- Com- 
parison of the first half (top panel) with the second half 
(bottom panel) indicates that the remapping process dif- 
fuses the discontinuity over a small number (~4) of zones 
initially, and "diffusion" slows substantially after the ini- 
tial phase, maintaining a somewhat consistent width in 
the discontinuity. 

9.6. Composition Remap 

Barring any nuclear or chemical transformations, the 
composition equation, eq. ([5]), states that Xi is con- 
served. The Lagrangian hydro portion of our algorithm 



will not change a cell's composition. All alterations in 
composition are, therefore, a result of advection, and in 
an ALE code advection is handled by the remapping al- 
gorithm. 

Equation ([5|) may be written in Eulcrian form: 



djpX, 
dt 



+ V • ipX.v) = , 



(144) 



which implies a close dependence of the advection of Xi 
on p and v. In practice, if the composition is not ad- 
vected in a manner entirely consistent with the advec- 
tion of mass, then the advection of composition will de- 
velop peculiarities. Since the remapping of density is 
truly a remapping of the subcell masses, we have de- 
signed a scheme for the remapping of composition which 
operates on an equivalent subcell quantity. 

Specifically, we define new Lagrangian quantities for 
composition. They are the number of species i in the cell, 
Ni^z, and the number of species i in the subcell, N[p. 
Analogous to the mass density is the number density, 
Ui^z and nfp. These new quantities are related to the 
previously defined quantities by: n,; — pXi and iV^ = UiV 
or Ni = mXi. 

The remap (or advection) of composition follows the 
scheme outlined for the advection of mass on the sub- 
cell level, with Ui replacing p and Ni replacing mass. 
The new compositions are then determined hy Xi^z — 
Ni^z/'nT'z- The final minor, but crucial, difference be- 
tween mass and composition remapping occurs during 
the repair process. While the repair process maintains 
the bounds on the number density, we also enforce bound 
preservation of the compositions. 

While we do not address specific rate equations in this 
work, transport and nuclear processes will require con- 
sideration of nonzero terms on the RHS of eqs. (O and 
(jl44p . Since our division of composition into subcell com- 
positions has consequences for discrete implementations 
with rates, we include here some discussion. 

Normally, the discrete form of eq. ([5]) with a nonzero 
RHS would be designed to operate at the cell level. In 
other words, it would involve Nza and Uza- However, the 
fundamental Lagrangian subunit is the subcell. There- 
fore, we have developed an algorithm to handle the com- 
position changes, AX^, at the subcell level due to changes 
at the cell level. 

Irrespective of the scheme employed, a condition which 
must be satisfied is that 



(145) 



z pGS(z) 



We'd like to convert the statement of number conserva- 
tion into the more useful relationship: 



E ^K,- 

pes{z) 



(146) 



When the change in composition due to the rates is ap- 
plied in operator-split form to the discrete hydro equa- 
tions, we have the following relations: 



and 



AiV, . = /^X,zmz , 



AiV^^p = ^Xl^m; , 



(147) 
(148) 
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Our dilemma is that the rate equations wiU determine 
AXz,i, but they put no constraint on AX^j. A fairly 
natural and simple choice is AX^ ^ — AX^^i, giving 

ATV;^, = AX,,,m^ . (149) 

Substituting eq. (|149p into cq. (|146p wc have 

(150) 

Hence, by simply stating that AX^^ — AXz^i, we have 
developed discrete rate equations that satisfy conserva- 
tion of species number and operates at the subcell level. 

9.7. Angular Velocity Remap 

The remap of the angular momentum deserves special 
mention. As in mass remapping, we try to find new sub- 
cell angular momenta, J^, based upon an approximate 
swept-region approach. In other words, 

eeS(z) 

which parallels eq. (|116p in form. A notable difference, 
however, is the expression for the swept angular momen- 
tum. It no longer involves a simple volume integral, but 
an integral weighted by (see appendix iBl for calculat- 
ing the discrete analog of this integral): 

SJe =^P j r^dV. (152) 

After finding the new subcell angular momentum we de- 
termine the new cell angular momentum, 

pes(z) 

and node angular momentum, 

Jp= E (154) 

zes(p} 

and in turn determine the new angular velocity; 

C^p = ^. (155) 

10. CODE TESTS 

In this section, we characterize BETHE-hydro's per- 
formance using several test problems. First, we demon- 
strate that this code produces 2"''-order accurate solu- 
tions for self-gravitating flows. To assess the accuracy of 
high Mach-number simulations, we use the Sod shock 
tube problem, Sedov blast wave, and Noh implosion 
problem, which all have analytic solutions. Furthermore, 
we simulate the Saltzman and Dukowicz piston prob- 
lems to test the code's ability to simulate piston-driven 
shocks using oblique meshes. In simulating two impor- 
tant hydrodynamic instabilities, the Rayleigh- Taylor and 
Kelvin-Helmholtz instabilities, we further demonstrate 
this code's strengths and limitations. Demonstrating 
BETHE-hydro's ability to simulate astrophysical phe- 
nomena, we conclude with a core-collapse supernova sim- 
ulation. 



10.1. ^'^-Order Accuracy 

To verify the 2"'i-order character of BETHE-hydro, a 
smooth hydrodynamics flow is required. The Goldreich- 
Weber self-similar collapse (see and Fig. [5]) satis- 
fies this requirement and is ideal to check convergence of 
the hydrodynamic and gravity solvers and their coupling. 
We use L^-nornis of the error: 

ii^^le.Ar'^l, (156) 

where the error is e^ = Pa,n&{xz) - Pz, Puns, is the ana- 
lytic density at t = 130 ms, Ar is the zone size, and d 
is 1 for ID and 2 for 2D. Strictly speaking, the simula- 
tions are Lagrangian and have time-varying zone sizes. 
Fortunately, this collapse problem is self-similar, imply- 
ing a direct correlation between the starting resolution 
and the resolution at a later time. Therefore, we use 
the initial zone size in eq. (|156p and in Fig. [T^ The 
L^-norm as a function of Ar (crosses) is plotted in Fig. 
[n[ As the figure demonstrates, both ID simulations (top 
panel) and 2D simulations (bottom panel) converge with 
roughly 2"'^-order accuracy. 

10.2. Sod Shock Tube Problem 

The Sod shock tube problem is a standard analytic 
test that assesses the code's ability to simulate three 
characteristic waves: a shock, a rarefaction wave, and a 
contact discontinuity. A simple gamma-law EOS is em- 
ployed with 7 = 1.4. Initially, the computational domain 
is divided into left and right static regions with different 
densities and pressures. Specifically, the left and right re- 
gions have densities of 1.0 g cm"'' and 0.125 g cm"'' and 
pressures of 1.0 erg cm"'' and 0.1 erg cm"^, respectively. 
This gives rise to a self-similar solution involving a shock 
propagating to the right, a rarefaction wave propagating 
to the left, and a contact discontinuity in between. 

In Fig. [T31 we display the Sod shock tube test re- 
sults for ID Lagrangian, ID Eulerian, 2D Lagrangian, 
and 2D Eulerian configurations. Cell-centered densities 
and locations are marked with plus signs, while the an- 
alytic results are denoted by solid dark lines. The ID 
Lagrangian result shows the appropriate density profile, 
while the other profiles have been shifted vertically so 
that distinguishing features are more easily compared. 
In addition, the profiles are further distinguished by dis- 
playing them at different times: ID Lagrangian (t — 0.2 
s), ID Eulerian (t = 0.225 s), 2D Lagrangian {t = 0.25 
s), and 2D Eulerian (t = 0.275 s). The ID calculations 
are resolved with 400 zones. Similarly, the 2D tests are 
resolved along the direction of shock propagation with 
400 zones, and they are resolved in the perpendicular 
direction by 10 zones. 

The overall features of the shock, post-shock material, 
contact discontinuity, and rarefaction wave are repro- 
duced. For all scenarios, the shock is resolved within 
a few zones. For the ID Eulerian and Lagrangian cal- 
culations the post-shock density is good to ~ 0.05% 
and ^ 0.01%, respectively. Upon viewing the profiles in 
greater detail, there is a noticeable departure from the 
analytic solution at the tail of the rarefaction wave where 
the density dips below the expected value. This error in 
density at the tail of the rarefaction wave is ^ 2.5% for 
the Lagrangian simulation and ^ 1% for the Eulerian 
simulation. Since the contact discontinuity moves with 
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the flow speed, both Lagrangian calculations resolve the 
contact discontinuity exactly from one zone to the next 
and maintain this resolution throughout the simulation. 
The Eulerian calculations, not surprisingly, distribute 
the contact discontinuity over several (~4) zones. 

10.3. Sedov Blast Wave 

A classic test, the Sedov blast wave problem provides 
a quantitative measure of a code's ability to simulate 
spherical explosions. Initial conditions are set so that the 
total energy of the Sedov blast is 0.244816 ergs, po = 1-0 
g cm~^, eo = 1 X 10"^*^ ergs g~^, and the gamma-law 
EOS has 7 = 5/3. We compare the simulations with the 
analytic result (solid lines) for the following scenarios: 
ID Lagrangian, ID Eulerian, 2D Lagrangian using the 
butterfly mesh, and 2D Lagrangian using the spiderweb 
mesh. 

Results of the Sedov explosion are plotted in Fig. [HI 
The top panel compares the analytic and numerical den- 
sity profiles, while the bottom panel shows the relative 
differences between the numerical solutions and the an- 
alytic profile. To easily distinguish different runs and 
details, we plot the ID Lagrangian calculation at i = 0.4 
s, the ID Eulerian at t = 0.53 s, the 2D Lagrangian using 
the butterfly mesh at t = 0.66 s, and the 2D Lagrangian 
using the spiderweb mesh at t = 0.80 s. The ID calcu- 
lations are resolved with 400 zones. The spiderweb test 
has a total of 12,381 zones with 200 radial zones and a 
maximum of 64 angular zones, while the butterfly mesh 
has a total of 35,000 zones with effectively 200 radial and 
200 angular zones. Qualitatively, all simulations repro- 
duce the overall structure and position of the shock and 
post-shock flows. 

There are some practical issues about formulating the 
Sedov runs that are worth mentioning. For the ID sim- 
ulations, we found that simply placing the initial explo- 
sion energy in a small number of inner zones was ad- 
equate. On the other hand, the 2D simulations using 
non-spherical grids required more care. Simply deposit- 
ing all the energy within a small number of zones near 
the center led to severe grid tangling and distortions. We 
remedied this by initiating all profiles with the analytic 
solution at t = 0.001 s. 

The bottom panel of Fig. [14] emphasizes the quan- 
titative accuracy of the post-shock solutions. Plotted 
are the relative errors of the density, (p — Pana)/Pmax, 
versus radius, where p is the simulated density profile, 
pana IS the analytic profile, and Pmax is the maximum 
density of the analytic profile. The ID Lagrangian sim- 
ulation gives the best results with a maximum relative 
error near the shock of ~ 2%. In comparison, the ID 
Eulerian test gives ^ 4% deviation near the shock. The 
2D Lagrangian simulation using the spiderweb mesh has 
a maximum deviation similar to the ID Lagrangian sim- 
ulation, even though the ID simulation has radial zones 
with half the zone size. The simulation using the spi- 
derweb mesh shows some slight departure from spherical 
symmetry. Specifically, for most of the post-shock re- 
gion, the peak-to-peak variation of density is less than 
~ 1%, and near the shock the peak-to-peak variation of 
density reaches ^ 3%. The 2D Lagrangian simulation us- 
ing the butterfly mesh has a maximum density deviation 
of ^ 4% near the shock and a deviation from spherical 
symmetry of similar magnitude. 



10.4. Noh Implosion Problem 

Initial conditions for the Noh Problem are uniform den- 
sity po = 1 g cm~'^, zero (or very small) internal energy 
£0 = ergs, a convergent velocity field with magnitude 
vt) — 1 cm s^^, and a gamma-law EOS with 7 = 5/3. 
Subsequent evolution produces a symmetric self-similar 
flow including supersonic accretion, an accretion shock, 
and stationary post-shock matter in which kinetic energy 
has been converted into internal energy. The analytic so- 
lution for this problem is 

[ |po(l-^) ,0,t;o} ifr>r, 

where v is the velocity magnitude, the shock position is 
rg = Ust, the shock velocity is Us = ^(7 ~ 1)^0, and d 
is 2 for 2D Cartesian, 3 for 2D simulations using cylin- 
drical coordinates, and 3 for ID spherically symmetric 
simulations. 

In Fig. [iSl we present the density profiles at i = 0.2 
s for ID, spherically symmetric simulations of the Noh 
problem using four resolutions. The initial grid spans 
the range x & {0 : 1} cm and is evenly divided into 200, 
400, 800, and 1600 zones. The solid line is the analytic 
solution for 7 = 5/3. Higher resolution simulations more 
accurately capture the shock position and the post-shock 
density profile. However, all resolutions depart signifi- 
cantly from the analytic solution near the center. This is 
"wall he ating" and is a common problem for Lagrangian 
schemes (lRidedl2000D . 

The four panels of Fig. [THj present the density profiles 
for 2D simulations. Three of the four panels show the 
results using Cartesian coordinates. The top-left panel 
shows the results using a Cartesian grid with 100x200 
zones, the top-right panel shows the results using the but- 
terfly mesh with 22,400 zones, and the lower-left panel 
shows the results using the spiderweb mesh. The fourth 
panel, lower-right, shows the results using 2D cylindrical 
coordinates and a Cartesian mesh with 100x200 zones. 
Obvious is the fact that the mesh employed has conse- 
quences for the solution. Both the top-left and bottom- 
right panels indicate that using the Cartesian mesh for 
this problem produces fairly smooth results, with some 
asymmetry (~7%) in the post-shock region. The lower- 
left panel shows that using the spiderweb mesh pro- 
duces perfectly symmetric solutions except near the cen- 
ter where the deviation from symmetry is as large as 
^--^25%. Finally, simulations using the butterfly mesh, 
top-right panel, shows a similarly mixed capacity to pre- 
serve symmetry. In this light, it is important to choose a 
grid that minimizes numerical artifacts for the problem 
at hand. A task easily accomplished with the use of ALE 
methods. 

10.5. Saltzman Piston Problem 

A standard te st for arbitrary grid codes, the Saltzman 
piston problem (|MargoHnlll988l ) addresses the ability of 
a code to simulate a simple piston-driven shock using a 
grid with mesh lines oblique to the shock normal. The 
top panel of Fig. [TTI shows the grid with 100 x 10 zones. 
For a grid with x Ny nodes within a domain defined 
by a; e {0 : Xmax} and y G {0 : j/max}, where 1.0 
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and Xraax - 

by 

x = (i- 1) 



0.1, the X positions of the nodes are given 



+ [Ny ~ j) sin 



(158) 

where i £ {I : N^} and j G {1 : Ny}. At the piston, the 
left wall is moving at a constant velocity, 1.0 cm s""'^, to 
the right. Initially, the density and internal energy are 
set equal to 1.0 g cm~'^ and 0.0 ergs, respectively, and 
we use a gamma-law EOS with 7 = 5/3. 

In the bottom panel of Fig. [T71 we show the grid and 
the density colormap at i = 0.925 s after the shock has 
traversed the domain twice, reflecting off the right and 
left walls once. In this snapshot, the shock is travel- 
ing to the right. Our results are to be compa red with 
Figs. 15 and 16 oflCampbell fc Shashkovl (|200ll) . Figure 
15 of lCampbell fc Shashkovl ! 2001h depicts a simulation 
with severe grid buckling, which we do not observe in 
our simulations. Instead, our Fig. 1 171 shows reduced grid 
buckling and appropria te densit ies in accordance with 
the results of Fig. 16 of lCampbell fc Shashkovl (|200lf ). 

10.6. Dukowicz Piston Problem 



Th e Dukowicz piston problem (Dukowicz fc Meltj 
|1992|) is another test using an oblique mesh. The ini- 
tial setup involves two regions. Region 1 has a density of 
1 g cm~'^ and is resolved with 144 x 120 zones. In the ver- 
tical domain, y € {0 : 1.5} cm, and the mesh lines evenly 
partition the space into 120 sections. Dividing region 1 
horizontally involves mesh lines with changing orienta- 
tion. The leftmost mesh line is parallel to the vertical, 
while the rightmost mesh line is oriented 60° relative to 
the vertical. The mesh lines in between smoothly transi- 
tion from 0° to 60°. Region 2 has a density of 1.5 g cm~'^ 
and is gridded with a 160 x 120 mesh, with the vertical 
mesh lines uniformly slanted at 60°. Including region 
1 and region 2, the bottom boundary spans the range 
X G {0 : 3} cm and is evenly divided into 304 segments. 
Initially, both regions are in equilibrium with P = 1.0 
erg cm~'^. The top, bottom, and right boundaries are 
reflecting, and the left boundary is a piston with a veloc- 
ity in the positive x direction and a magnitude of 1.48 
cm s~^ (see the top panel of Fig. |TS]for a low resolution 
example of this grid). 

A piston-driven shock travels from left to right, and 
encounters the density jump at an angle of 60°, produc- 
ing a rich set of phenomena. The incident shock contin- 
ues into the lower density region, a transmitted/refracted 
shock propagates into the higher density region, a vor- 
tex sheet develops behind the transmitted shock, and a 
reflected shock propagates into the incident shock's post- 
shock flow (see the labels in the bottom plot of Fig. [TS] 
for visual reference). The results are to be compared 
with th e semi-analytic shock-polar analysis presented by 



Dukowi cz fc Mcltz (199^ (see Figure 13 and Table I of 
Dukowicz fc Meltd (|1992l) ). In their paper, angles sub- 



tended by the five regions are presented. In Table |21 we 
recast this information as the angles that the transmitted 
shock, vortex sheet, reflected shock, and incident shock 
have with a:-axis. The first row gives the analytic values, 
which should be compared with the simulated orienta- 
tions in the second row. Strikingly, despite the fact that 
the incident shock has traversed a grid with an oblique 



mesh, the simulated and analytic orientations differ very 
little. The simulated and analytic reflected-shock orien- 
tations agree to within ~ 0.5%, and the orientations of 
the transmitted shock and vortex sheet agree to within 
2%. 

10.7. Rayleigh- Taylor Instability 

Here, we simulate the growth of a single mode of the 
Rayleigh- Taylor Instability. A heavy fluid is placed on 
top of a lighter fluid in the presence of a constant gravi- 
tational acceleration. The top and bottom densities are 
pi = 2.0 g cm~'^ and p2 = 1.0 g cm"'^, respectively, 
and the gravitational acceleration points downward with 
magnitude g = 0.1 cm s~^. The domain is a rectangle 
with X e {-0.25 : 0.24} cm and y € {-0.75 : 0.75} cm 
and 100 X 300 grid zones. The top and bottom boundaries 
are reflecting while the left and right boundaries are peri- 
odic. For such configurations, small perturbations of the 
interface between the heavy and light fluids are unstable 
to exponential growth. Assuming that the boundaries 
are far from the interface, the exponential growth rate of 
a perturbation with wavenumber k is 



I kgjpi - P2) 
pi + P2 



(159) 



For all single-mode Rayleigh- Taylor simulations we initi- 
ate the perturbation by setting the vertical component of 
the velocity equal to Vy = 2.5 x 10~^(l-|-cos(27ra;/A))(l-|- 
cos(37r2/)), where the wavelength. A, is 0.5 cm. Therefore, 
we simulate one wavelength of the mode, and the expo- 
nential growth rate (w) should be 0.65 s~^. 

Figure [19] shows the evolution of a single-mode 
Rayleigh- Taylor instability dX t — 12.75 s for four res- 
olutions. From left to right, the grid sizes are 50 x 150, 
74 X 222, 100 X 300, and 150 x 450. Gross features 
compare well, with all the resolutions differing by only 
a few percent in the maximum and minimum position 
of the interface. However, higher resolution simulations 
manifest greater complexity for the Kelvin-Helmholtz 
rolls. One can compare the thir d pane l to the results 
of Fig. 4.5 of iLiska fc Wendrofj (|2003l) . In a general 
sense, these authors conclude that fewer features imply 
more dissipation. However, it could be that some of the 
Kelvin-Helmholtz rolls are seeded by grid noise in some 
schemes^. 

Next, we discuss the effects of artificial viscosity on 
the single-mode Rayleigh- Taylor flow. As stated in 
artificial viscosity is a requirement for ALE algorithms 
in order to simulate shocks. In the current formula- 
tion, there are two parameters of the artificial viscosity 
scheme. One parameter, c2, is the coefficient for (V • z/)^, 
which is largest in shocks. Effectively, this term provides 
shock resolution. The second parameter, ci, is the coef- 
ficient of the term that is proportional to Cs(V • v) and 
is designed to reduce the amount of post-shock ringing 
in t he solution. Typical v alues suggested for both are 
1.0 (|Campbell fc Shashkovl lIOOl]. For the single-mode 
Rayleigh- Taylor test, the flows are subsonic, so the term 
multiplied by ci has the greatest impact on the mag- 
nitude of the artificial viscosity forces. Each panel of 

^ See results from Jim Stone's Athena for a more fa- 
vorable comparison and a discussion of the grid nois e issue 
l |http://www.astro.princeton.edu~jstone/tests/rt/rt.html1 . 
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Fig. [20] presents the nonhnear Rayleigh- Taylor flow af- 
ter 12.75 s for different values of the artificial viscos- 
ity parameter, ci. Clearly, ci = 0.01 (left panel) re- 
produces the expected results for this low-Mach-number 
flow. On the other hand, the ci = 0.1 run (center panel) 
displays significant departures for the Kelvin-Helmholtz 
rolls, while the overall progression of the plumes remains 
similar. Unfortunately, the model with ci = 1.0 (right 
panel) completely suppresses the Kelvin-Helmholtz rolls 
on these scales and plume progression is severely re- 
tarded. 

In Fig. [211 the interface perturbation amplitude is 
plotted as a function of time (dashed lines). The solid 
line is the analytic exponential growth rate scaled to the 
simulation results. The dashed lines show the simulation 
results for viscosity parameters, ci, of 0.01, 0.1, and 1.0. 
Three distinct phases are apparent: an early transient 
phase, a phase in which the slope most closely matches 
the exponential growth rate, and the subsequent non- 
linear phase. The simulations with ci — 0.01 and 0.1 
manifest exponential growth for several e-folding times. 
The run with ci = 1.0, on the other hand, seems to fol- 
low the linear phase for only 1 s (~ 1/2 e-folding), if at 
all. Around 5 s, the evolution of the interface amplitude 
enters the nonlinear phase. 

For the artificial viscosity scheme that we employ, the 
Rayleigh- Taylor instability test indicates that lower val- 
ues of ci are preferred. However, tests of the Sod shock- 
tube problem including a parameter study of ci indi- 
cate that unwanted post-shock ringing is diminished only 
for values above ^ 0.5. This represents the primary 
weakness of BETHE-hydro to simulate fiows with both 
shocks and Rayleigh- Taylor instabilities. To be clear, 
this does not represent a weakness of ALE methods in 
general, but of the tensor artificial viscosity algorithm de- 
signed to mitigate the hourglass instability that we em- 
ploy. Othe r artificial viscosity sch e mes s uch as edge vis- 
cosity (see iCampbell fc Shashkovl (|2001[ ) for references) 
allow for proper development of hydrodynamic instabil- 
ities, but do little to suppress the hourglass instability 
(Milan Kucharik, private communication). Of great in- 
terest to users of ALE is a methodology that suppresses 
the hourglass instability and post-shock ringing while en- 
abling proper evolution of hydrodynamic instabilities. 

10.8. Kelvin-Helmholtz Instability 

Another important phenomenon we explore is the 
Kelvin-Helmholtz shear instability. lAgertz et al.l (|2007l ) 
have shown that SPH has trouble simulating the Kelvin- 
Helmholtz instability when extreme density contrasts are 
involved. We find that this instability is reasonably well 
handled in BETHE-hydro, and that the evolution during 
the small amplitude regime is accurately characterized by 
analytic linear analysis. The calculational domain cov- 
ers the square region bounded by a; € {0 : 1} cm and 
y G {0 : 1} cm and has 256 x 256 zones. The top and 
bottom boundaries are reflecting, while the left and right 
boundaries are periodic. For y < 0.5 cm, pb = 1-0 g 
cm""^, and for y > 0.5 cm, pt — Pb/x, where x — ^■ 
These regions are in pressure equilibrium with P = 1.0 
erg cm~^. The shearing velocity, Wshear is scaled with 
respect to the slowest sound speed (sound speed in the 
bottom region, Cfc). This relative shearing velocity is split 
between the top region, which flows to the left with speed 



(160) 



ji'shcar, and the bottom region, which flows to the right 
with the same speed. For this test, the linear coefficient 
in the viscosity, ci, is set to 0.01. A gamma-law EOS is 
used with 7 = 5/3. The initial perturbation (see Agertz 
et al. 2007 ) is placed in a small band centered on the 
interface and is a perturbation in velocity given by 

/ 2ttx\ 

Vy = SvyVshcar siu ( j for | J/ - 0.5 1 < 0.025 , 

where A = 1/3 cm. It should be noted that this pertur- 
bation is not an eigenmode of the instability. Therefore, 
simulations have a transient phase at the beginning in 
which this perturbation settles into one or more modes 
of the instab ility. We have found that using amplitudes 
suggested bv lAgertz et~all (|2007D . Svy = 1/80 and 1/40, 
produces strong transients that complicate the interpre- 
tation of the linear regime. To avoid strong long-lasting 
transients, we set dvy = 1/160. These initial conditions 
should produce a perturbation of the interface that grows 
exponentially in magnitude with an e-folding time given 

by 

'^(Ptop + Pbot) (161) 



Tkh 



27r'y shear -/PtopPbo 



For the simulations presented here, we consider two 
shearing velocities: Wshcar — \cb = 0.6455 cm s~^ and 
Wshear = = 0.32275 Cm s~^. Hence, tkh is 0.262 s 
and 0.523 s, respectively. 

The top panel of Fig. [^ shows the evolution at t = 5.5 
s for t^shear = \cb- The first set of nonlinear Kelvin- 
Helmholtz rolls have appeared. Qualitatively, the mor- 
phology of the rqlls is similar to the results presented in 
Figure 13. of lAgertz et all labeled by "grid 1" at 

t — 27rrKH (there is a factor of 2tt difference in the def- 
inition of TKH between our work and theirs). There are 
two main differences. For one, we simulate with a factor 
two larger wavelength to adequately resolve the wave- 
length for linear analysis. Secondly, the time at which 
we present the results in the top panel of Fig. [221 cor- 
responds io t — 1.67 X (27rrKH), not t = 27rrKH- We 
have noticed similar discrepancies in the time required to 
achieve similar evolution using A = 1/6 and Ushoar = \cb. 
Additionally, we analyze the growth rate of the interface 
during the linear regime. We determine the interface 
amplitude by generating a contour for p = 0.5(pt + Pb), 
measuring the peak to peak amplitude, and dividing by 
two. The bottom plot of Fig. [22| shows this amplitude 
(solid line) versus time and compares to the expected 
exponential growth (dashed line) for a simulation with 
"shear = \cb (green) and Ushear = \cb (blue). There 
are three distinct phases in the log- linear plot: an early 
transient phase, a phase in which the slope most closely 
matches the exponential growth rate, and the subsequent 
nonlinear phase. While the growth rate during the linear 
regime differs slightly (~10%) from theory, the bottom 
panel of Fig. [221 demonstrates that the correct depen- 
dence of the growth rate on Wshear is obtained. 

10.9. Core- Collapse Test 

Incorporating all aspects of BETHE-hydro, we sim- 
ulate ID and 2D hydrodynamic core-collapse of a 
15 Mp, star. We us e th e sl5s7b2 model (S15) 
of IWooslev fc Weaved (flQQSh and the Shen EOS 
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(|Shen et al.1 |1998[ ) and initiate collapse usi ng a Y^- 
p parameterization (jLiebendorfer et al] l2006f ). While 
there are no analytic solutions for such a test, we 
test to see whether the density profiles, timescales, 
and shock radii, etc. all matc h experience with 
other cod es and results p u blishe d (iLiebendorfer et aTl 
2001bllal: [Rampp k Jankal 120021: iBuras et all l2003l: 
Thompson et al.ll2003l: ILiebendorfer et al.ll2005D For 2D 
simulations, the grid is composed of a butterfly mesh in 
the interior with a minimum cell size of ~0.5 km and 
extends to 50 km where a spherical grid carries the do- 
main out to 4000 km. In total, there are 23,750 cells, 
with an effective resolution of ^^250 radial and ~100 an- 
gular zones. For the best comparison, the ID grid has 
250 zones, mimicking the effective radial resolution of the 
2D simulation. Figure [531 depicts the density vs. radius 
for ID (lines) and 2D (crosses) at times 0, 70, 110, 130, 
140, and 150 ms after the start of the calculation. Core 
bounce occurs at 148 ms. Other than a ~10% difference 
in the shock radii at 150 ms, ID and 2D calculations 
track one another quite well. Furthermore, Fig. l24l shows 
an exceptional correspondence between the ID and 2D 
gravitational accelerations. 

Using this core-collapse test to best represent the 
conditions during astrophysical simulations, we describe 
some timing results of this problem. A standard mea- 
sure is the average CPU time spent per cell per cy- 
cle. Using one core of a Dual-Core AMD Opteron™2.8 
GHz processor for this test, the average time spent is 
8.125 X 10"^ CPU seconds/cycle/cell. 80% of which is 
spent in AMG1R6, the multigrid solver. For this prob- 
lem, the iterative multigrid solver usually takes ~5 to 
'^lO cycles to achieve a fractional residual of ^--^10"^ to 
~10~^, where the fractional residual of the linear system 
Ax = h is {Ax - b) /b. 

We re-simulate the collapse of the S15 model, but 
with an angular velocity profile of the form Q{r) = 
1/(1 + (r/A)2), where A = 1000 km and no = 2 ra- 
dians s~^. The vs. r plot is shown in Fig. [25l Core 
bounce occurs at 151 ms. First of all, total angular mo- 
mentum is conserved to machine accuracy. Since there 
are no analytic descriptions for the evolution of angular 
velocity in core collapse scenarios, we cannot validate the 
numerical results via analytic theory. However, the cen- 

2/3 

tral angular velocity, flc is proportional to pc , where 
Pc is the central density. Since the central density com- 
presses from ^10^° g cm^"^ at t = ms to ^2.2 xlO^"* g 
cm~^ at t = 160 ms, ft^ should be ^1600 radians s""'^ at 
t ~ 160 ms. This is consistent with results shown in Fig. 
1251 Furthermore, the angular velocity evolves smoothly 
with no evidence of axis effects. On the other hand, there 
is a slight, but noticeable, glitch at the location of the 
shock. 

11. DISCUSSION AND CONCLUSIONS 

In this paper, we have presented the algorithms em- 
ployed in BETHE-hydro, a new code for ID and 2D 
astrophysical hydrodynamic simulations. The hydrody- 
namic core is an ALE algorithm, and its most strik- 
ing feature is the ability to use arbitrary, unstruc- 
tured grids. With finite- differencing based upo n the 
su pport-operator in e thod (iShashkov fc Steinb erg 1995) 
of ICaramana et aTl (|1998D and ICaramana fc Shashkov. 
(|1998f) . energv is conserved to roundoff error in the ab- 



sence of rotation and gravity, and momentum is strictly 
conserved using Cartesian coordinates. For all other cir- 
cumstances, energy and momentum are conserved ac- 
curately, if not precisely. We use a subcell remapping 
scheme that conservatively remaps mass, momentum, 
energy, and number density. For 2D calculations us- 
ing cylindrical coordinates, we include rotational terms 
in the Lagrangian solver, and develop a remapping al- 
gorithm that conservatively remaps angular momentum 
while minimizing unwanted features in the angular ve- 
locity near the axis. To provide shock resolution and 
grid stability, we use t he te nsor artificial viscosity of 
iCampbell fc Shashkov! ([200lh . and to minimize hour- 
glass instabilities, we have developed a subcell pressure 
me thod that is a close deriv ative of the scheme developed 
by ICaramana fc Shashkov! ([1998) , but avoids pathologi- 
cal problems when used in conjunction with the subcell 
remapping algorithm. Finally, we have developed a grav- 
ity solver for arbitrary grids that uses a support-ope rator 
technique for elliptic equations (iMorel et al.lll998l ) and 
an iterat ive multigrid-precond itioned conjugate-gradient 
method (jRuge fc Stuben|[l987D to solve the system of fin- 
ear equations. 

Overall, BETHE-hydro offers many unique and useful 
features for astrophysical simulations. For one, by us- 
ing ALE techniques, the structure of BETHE-hydro is 
straightforward, enabling simple inclusion of a variety of 
additional physics packages. Examples, which will be dis- 
cussed in future papers, are nuclear networks and time- 
dependent radiation transport. In contrast with other 
techniques, such as higher-order Godunov methods, no 
assumptions are made in ALE techniques about char- 
acteristic waves nor the relationships among important 
thermodynamic variables. Hence, one of its useful fea- 
tures is an ability to use a general EOS. 

Most important among BETHE-hydro's strengths is 
the ability to solve self-gravitating hydrodynamic flows 
on arbitrary grids. This is achieved primarily because 
the foundation of the Lagrangian hydrodynamics solver 
is an arbitrary, unstructured polygonal grid. Further- 
more, grids may be time-dependent since the hydrody- 
namic flow from one timestep is remapped to another 
arbitrary grid for subsequent evolution. Consequently, 
simulations may be executed using a purely Lagrangian, 
purely Eulcrian, or an arbitrarily defined time-dependent 
grid. With BETHE-hydro, simulations have great flexi- 
bility, tailoring the grid to minimize numerical error and 
suiting the grid to the computational challenge. 

Ironically, this flexibility leads to BETHE-hydro's most 
prominent weakness, which is shared among all ALE 
codes, the hourglass instability and grid buckling. We 
use subcell pressure and tensor artificial viscosity algo- 
rithms to mitigate these numerical artifacts, but the res- 
olution is imperfect. In SjHl we present a subcell pressure 
algorithm that is compatible with the subcell remapping 
algorithm, but the efficacy of hourglass elimination is 
slightly compromised. Although the tensor artificial vis- 
cosity does well to mitigate grid buckling, it can be more 
resistive to the proper development of hydrodynamic in- 
stabilities compared to other artificial viscosity schemes 
used in ALE (Milan Kucharik, private communication). 
Hence, we diminish the effects of the hourglass instabil- 
ity, but with some unwanted side effects. While mitiga- 
tion in 2D is tractable, the instability in 3D has many 
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more modes and is not as easily eliminated (Guglielmo 
Seovazzi, private communication ), making long 3D ALE 
simulations a challenge. 

Throughout this paper, we have not only demonstrated 
BETHE-hydro's flexibility, but have shown that it pro- 
duces accurate and 2"'^-order convergent solutions. With 
density distributions having analytic potentials, we have 
shown that the 2D gravity solver gives accurate spheri- 
cally symmetric potentials when a non-spherical grid is 
used, it produces accurate non-spherical potentials, and 
solutions converge with 2"'^-order accuracy. Further tests 
demonstrated accurate solutions for self-gravitating hy- 
drostatic and dynamic problems. We have shown that 
our hourglass elimination algorithm minimizes the hour- 
glass instability and does not present problems when 
used in conjunction with subcell remapping. In addition, 
wc have quantified the accuracy of hydrodynamic simu- 
lations by simulating problems with analytic solutions. 
Simulating piston-driven shocks using oblique meshes, 
we confirmed we can obtain accurate solutions in the 
context of arbitrary grids. Verifying the code's ability 
to capture basic hydrodynamic instabilities, wc simulate 
the Rayleigh- Taylor and Kelvin-Helmholtz instabilities. 
Concluding the tests, we simulated a supernova core col- 
lapse, which demonstrates the ability to simulate com- 
plex astrophysical phenomena. 

Simulating hydrodynamic flow is fundamental to un- 



derstanding most astrophysical objects, and despite the 
long tradition of hydrodynamic simulations many puz- 
zles remain. This is due primarily to the need to ad- 
dress time-dependent gravitational potentials, compli- 
cated equations of state (EOSs), flexible grids, multi-D 
shock structures, and chaotic and turbulent flows. There- 
fore, with BETHE-hydro, we introduce a uniquely flexi- 
ble and functional tool for advancing the theory of com- 
plex astrophysical phenomena. 



We would like to thank Mikhail Shashkov, Raphael 
Loubere, Milan Kucharik, and Burton Wendroff for 
valuable conversations about ALE techniques. Dis- 
cussions concerning the 2D gravity solver with Ivan 
Hubcny, Jim Morel. David Moulton, and David Keycs 
of TOPS were very fruitful. We acknowledge support 
for this work from the Scientific Discovery through Ad- 
vanced Computing (SciDAC) program of the DOE, un- 
der grant numbers DE-FC02-01ER41184 and DE-FC02- 
06ER41452, and from the NSF imdcr grant number AST- 
0504947. J.W.M. thanks the Joint Institute for Nu- 
clear Astrophysics (JINA) for support under NSF grant 
PHY0216783.' We thank Jeff Fookson and Neal Lauber 
of the Steward Computer Support Group for their in- 
valuable help with the local Beowulf cluster, Grendel. 



APPENDIX 

HYDRODYNAMIC EQUATIONS: EULERIAN FORM 
In Eulerian form, the equations for conservation of mass, momentum, and energy are 

dp 



dt 



+ V-(pv) = 0, 



and 



9{pv) 

dt 

di 



+ V • {pvv) -h VP = -pV$ , 



-I- V • [{pe + P+ l/2pv^)v\ = -pV$ • V. 



(Al) 
(A2) 

(A3) 



VOLUME AND WEIGHTED- VOLUME INTEGRALS IN 2D 



In this appendix, we give the exact analytic volumes and weighted-volume integrals in discrete form. Since we are 
dealing with 2D simulations all integrals are 2D integrals, and can be further reduced using Green's formula to ID 
boundary integrals. For example, the volume of zone z in Cartesian coordinates is 



or alternatively, 



/ dV = (f xdy = ^l-{xi+X2){y2-yi) 



/ dV = -i ydx = -Y\\{yi + y2){x2 - xi) , 



(Bl) 



(B2) 



where the direction of the boundary integral is clockwise, e indicates an edge of the cell, and subscripts 1 and 2 
represent the endpoints of edge e. In cylindrical coordinates, the volume integral is 



[ dV = 27rJ / r'^dz = 27r V i(r? -h nrs + rl){z2 - z{) . 



(B3) 



Throughout BETHE-hydro, integrals of linear functions produce volume integrals weighted by x or y. The discrete 
form of these integrals are 



/ xdV=\i x^ dy = ^S\\{x\ ^ xxX2 ^ x\){y2 - y\) 



(B4) 



25 



and 



/ ydV = -l(f 2/^rfa; = -Vi(yf + 2/12/2 + yi)(a;2-a;i). 



Using cylindrical coordinates these are 

/ rdV = 2itI (f r^dz = 27r V -^(r? + rfrs + nr^ + 4){z2 - zi) 
Jv^ ^ Jav^ g ^■^ 

and 

Jy zdV = 27r^ r'^zdz 

= 277 Eeg[(H + 5^2ri + H)(^2-^i) + (ri+r2ri+rf)zi] (2:2-2:1). 

Finally, for angular momentum remapping, we require integrals in cylindrical coordinates weighted by r^: 
/ r'^dV = 2tt- <f> r"^ dz = 2t:S^ ^{rf + rlr2 + rlrl + rxrl + 4){Z2 - zi). 

JVz 4 Jqy^ ^ M 



(B5) 

(B6) 
(B7) 

(B8) 
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TABLE 1 

The predictor-corrector procedure. The operations to determine 
predicted and corrected valu es ar e listed in the left and right 
columns, respectively. see ^4.21 for, a detailed discussion. 



Predictor Corrector 



pn ^ gn > pi pn+l/2 ^ _gn + l/2 ^ ^1 + 1/2 

/'I'isc. r+'''^ /^isc fer+'Z' — * '^+' 

^ + 1/2 _ 1 J-^+l.pr _j_ ^1 + 1/2 _ 1 (-^+1 -(_ 

/*" ^{j« + l/2 > ^n + l.pr jn + 1/2 ^^ + 1/2 ^ ^n+l 

^n + l.pr ^ ^n + l,pr ^ ^n + l.pr ^-" + 1 ^ ^^ + 1 ^ p'^^"'- 

calculate & A^+^'f 

gn + l,pr^pn + l,pr ^ ^n4-l,pr 

+ _ l^^n + l,pr _|_ 

j3n-|-l/2 _ 1 ^pn-\-l,pr _|_ ^tiJ 

g" + l/2 = l(g« + l,pr 

7-71 + 1/2 _ 1^ 7>n + l,pr I Xn\ 
^p — 2 l^P + ^p J 



TABLE 2 
Relative error in density, 

(pz - Pana)/(max(pana)), AT THE CENTER 

FOR THE GOLDREiCH- Weber self-similar 

COLLAPSE TEST. THE FIRST COLUMN IS THE 
TIME IN MILLISECONDS, AND THE SECOND 

AND THIRD COLUMNS ARE THE ERROR 
BETWEEN THE SIMULATION, Pz, AND THE 
ANALYTIC SOLUTION, Pana, SCALED BY THE 

MAXIMUM DENSITY OF THE ANALYTIC 
SOLUTION max(pana) FOR ID SIMULATION 
AND 2D SIMUL ATIONS, RESPECTIVELY. SEE 
^5.5| F0R A DISCUSSION. 



Time (ms) 


ID 


2D 







-9.2 X 10-** 


-1.2 X 10" 




20 


-7.2 X 10-"* 


-1.4 X 10" 


-3 


40 


-1.2 X 10-3 


-2.2 X 10" 


-3 


60 


-1.9 X 10-3 


-3.6 X 10- 


-3 


80 


-3.2 X 10-3 


-6.0 X 10" 


-3 


100 


-6.2 X 10-3 


-1.1 X 10" 


-2 


120 


-1.6 X 10-2 


-2.8 X 10" 


-2 


130 


-3.6 X 10-2 


-6.3 X 10" 


.2 



TABLE 3 

The analytic and simulated orientatio ns of features for the Dukowicz 
PROBLEM. While [Dukowicz fc Melt3 II1992I) published the angles subtended by 

VARIOUS REGIONS, WE LIST THE ANGLES OF EACH FEATURE WITH RESPECT TO THE X-AXIS. 



Transmitted Shock Vortex Sheet Reflected Shock Incident shock 



Analytic (deg.) -101.02 -122.92 146.12 90 

Simulation -103.21 -124.42 146.92 90.32 
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Fig. 1. — The flowchart for BETHE-hydro. See i|2]for a discussion. 
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Fig. 2. — On the left, a representative butterfly mesh and on the right, a representative spiderweb mesh. In general, we eonstruct meshes 
with much higher resolution. However, to easily see the features of the mesh we reproduce low resolution versions here. The butterfly 
mesh is a standard grid that captures the benefits of a spherical grid near the outer boundary and avoids the singularity at the center. 
The spiderweb mesh captures the benefits of a spherical grid throughout, but near the center, the angular resolution is modified to avoid 
extreme Courant-condition limits and the singularity. Construction begins with half of an octagon at the center. Subsequent tiers of nodes 
are placed Ar farther in radius from the interior set. When rAO exceeds Ar the angular size is halved. For these transition nodes, Ar is 
alternately multiplied by 1 + e and 1 — e to exaggerate the concavity of the cells. 
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Fig. 3. — A portion of an unstructured, polygonal mesh. Filled circles indicate node positions, p. Crosses mark cell centers, z. Open 
circles show the mid-edge locations. Cells are bounded by the edges (solid lines), and the cell is further divided into subcells via dashed 
lines. C'p-i- and Cp„ are the half-edge area vectors, and ap_j are the area vectors for the lines connecting the cell centers and mid-edges 
(see 33] for a discussion). 
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Fig. 4. — Plots of ti0 vs. r (top panel) and Q vs. r (bottom panel) for the simple rotational test discussed in ^4.3.21 We employ three 
remapping regions. The inner 0.1 cm is Eulerian, zones exterior to 0.2 cm follow Lagrangian dynamics, and the region in between provides 
a smooth transition between the two domains, and Q profiles are presented at t = 0.0, 0.1, 0.2, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, and 1.0 s. 
In the top panel, the deviation of the simulation (crosses) from the analytic solution (solid lines) is not discernible. The maximum error as 
measured by (j^^^ana ~ '''0)/niax(ti0 ana) is ^ 4 X 10~^ near the beginning of the simulation and ~ 8 X 10~* at the end. The C vs. r panel 
(bottom) shows similar accuracy, except at the center, where (Qana — f^) reaches a maximum value of 0.2 rads s^-"^ at t = 0.55 s. 
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Fig. 5. — Potential of a uniform density sphere. The density of the sphere is 1 g cm~ and the radius is 1 cm. Presented are results for 
four resolutions of the butterfly mesh (see left panel of Fig. [JJ: 2550 cells with an effective radial resolution Ar ~ 0.02 cm (blue); 8750 
cells with Ar ~ 0.01 cm (green); 15,200 cells with Ar ~ 0.006 cm (yellow); and 35,000 cells with Ar ~ 0.005 cm (red). The top panel 
shows (<& — 'I'ana)/|'I'ana|, where # are the cell-center potentials as determined by the Poisson solver, and #ana is the analytic potential. 
The bottom panel shows the minimum error for each resolution as a function of the effective radial resolution. The solid line illustrates the 
fact that the 2D Poisson solver converges with 2"'^-order accuracy. 
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Fig. 6. — The potential of a homogeneous spheroid with e = 0.8 (top panel) and the relative accuracy of the potential, — <I'ana)/|'3?ana| 
(bottom panel), po = 1 g cm~^ and the equatorial radius is = 1 cm. For the given eccentricity, the polar-axis radius, rt, is 0.6 cm. 
The grid is a butterfly mesh, but the outer boundary follows the ellipse defining the surface of the spheroid. With N^^n = 35, 000 and 
Ar/ra ^ 0.005, the relative error in the potential ranges from ~2 xl0~® near the outer boundary to ~3 xlO~^ in interior regions. Features 
in the relative error that track abrupt grid orientation changes in the mesh are apparent, but the magnitude of these features does not 
dominate the errors. 
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Fig. 7. — Hydrostatic Equilibrium. The grid is the butterfly mesh with 8750 zones. The initial model is a Lane-Emden polytropc with 
7 = 5/3, M = 1 Mq, and R = 2.9 X 10^" cm. The green crosses show the density profile at t = 1 X 10* s. The solid line shows the maximum 
density as a function of time. The oscillations are due to the slight difference between an analytic hydrostatic equilibrium structure and a 
discretized hydrostatic equilibrium structure. This figure shows the code's ability to follow oscillations for many periods with very little or 
no attenuation. 
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Fig. 8. — Density profiles for Goldreich- Weber self-similar collapse in ID (top panel) and 2D (bottom panel) and at t =0, 20, 40, 60, 
80, 100, 120, and 130 ms. Simulation results (crosses) are compared with analytic solutions (solid lines). A gamma-law EOS is used with 
7 = 4/3. The physical dimensions have been scaled so that M = 1.3 Mq, the initial central density is 10^" g cm~'^, and the maximum 
radius of the profile is 1.66 X 10* cm. The initial grid for the ID simulation uses 200 evenly-spaced zones. For the 2D simulation, a butte rfly 
mesh with 35,000 zones with effectively 200 radial zones is used. These tests are calculated using the Lagrangian configuration. See i|5.5l 
and Table [2] for quantitative discussion of the accuracy. 
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Fig. 9. — Comparison when hourglass suppression is turned off (left) and on (right) for the Goldreich- Weber self-similar collapse. The 
simulations in this figure were run in the Lagrangian configuration with the butterfly mesh. Displayed in the region "r > 0" are results at 
t = 0.118 s when the hourglass suppression scheme of ^is used and the scale factor has been set to 2.0. The results shown in the region 
"r < 0" represent what happens when the hourglass suppression scheme is turned off. 
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Fig. 10. — Single-mode Rayleigh- Taylor instability: investigating hourglass suppression schemes (see i|10.7l for a description of the setup). 
We prese nt the results of the single-mode Rayleigh- Taylor instability with no hourglass fix (left panel), the hourglass subcell pressure 
scheme of lCaramana &: Shas hkov (1998) (center panel), and the modified subpressure scheme described in this paper (right panel). While 
the results of the left panel arc acceptable, there are hints of hourglass patterns in the co ntours. The central panel is a consequence of the 
incompatible nature of the subcell remapping scheme and the subcell pressure scheme of lCaramana fc Shashkovl (119981 ). and our modified 
subpressure scheme (right panel) suppresses the hourglass distortions, while preserving the expected flow. 
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Fig. 11. Remapping a ID step function. A step function in density is remapped many times with a grid that oscillates for two full 
cycles. There are 51 nodes (50 cells) and rimax = 800 remapping steps. The top panel displays the density profile for remapping steps to 
l/2nrrnmax, and the bottom panel shows the profile for steps l/2nniax to rimax- Initially, the discontinuity spreads over a small number of 
zones (~4), and in later steps (bottom panel) the discontinuity spreads very slowly. 
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ir [cm] 

Fig. 12. — The L^-norm as a function of Ar (crosses), the initial zone size, is plotted for Goldreich- Weber simulations in ID (top panel) 
and 2D (bottom panel). The L^-norm is calculated at t = 130 ms for both simulations. Both ID and 2D simulations (crosses) converge 
with roughly 2"'^-order accuracy (solid line). Sec illO.ll 
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Fig. 13. — Sod shock tube test. For a gamma-law EOS and 7 = 1.4, we compare the results of the ID Lagrangian (bottom), ID Eulerian 
(2nd from the bottom), 2D Lagrangian (3rd), and 2D Eulerian (top) Sod shock tube tests with the analytic result (solid lines). Other than 
the ID Lagrangian results, the density profiles have been shifted vertically to distinguish features. The profiles are further separated by 
displaying them at different times: ID Lagrangian {t = 0.2 s), ID Eulerian {t = 0.225 s), 2D Lagrangian {t = 0.25 s), and 2D Eulerian 
{t = 0.275 s). The ID calculations are resolved with 400 zones, and the 2D tests are resolved with 400x10 zones. See i|10.2l for a discussion. 
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Fig. 14. — Density profiles of the Sedov blast wave test. The Sedov blast energy is 0.244816 ergs. Initial conditions are po = 1.0 g 
cm~^, £o = 1 X 10"^'' ergs g~^, and the gamma-law EOS has 7 = 5/3. Top panel compares the simulations (crosses) with the analytic 
result (solid lines) for ID Lagrangian [t = 0.4 s), ID Eulerian (i = 0.53 s), 2D Lagrangian using the butterfiy mesh (i = 0.66 s), and 2D 
Lagrangian using the spiderweb mesh [t = 0.80 s). ID calculations are resolved with 400 zones. The spiderweb test has a total of 12,381 
zones with 200 radial zones and a maximum of 64 angular zones. The butterfiy mesh has a total of 35,000 zones with effectively 200 radial 
and 200 angular zones. Plotted in the bottom panel are the relative errors of the density, (p — pana)/pmax, vs. radius, where p is the 
simulated density profile, pana is the analytic profile, and pmax is the maximum density of the analytic profile. 
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Fig. 15. — Density profile for tfie Noli Problem at t = 0.2 s. The solid line is the analytic solution for 7 = 5/3, and the crosses show the 
numerical results for resolutions of 200, 400, 800, and 1600 zones. The downturn of the density profile near the center is "wall heating," 
a common problem for Lagrangian schemes (Rider 2000). For the post shock material, the higher resolution runs capture the analytic 
solution. The upstream flow matches the analytic solution for all resolutions. 
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Fig. 16. — Similar to Fig. 1151 except here we present the results of 2D tests. All but the lower-right show results using 2D Cartesian 
coordinates. The lower-right shows results using 2D cylindrical coordinates. The grids used are: a Cartesian grid with 100x200 zones 
(top-left and lower-right), a butterfly mesh with 22,400 zones (top-right), and a spiderweb mesh with 8550 zones (lower-left). Both the 
top-left and bottom-right panels indicate that using the Cartesian mesh produces fairly smooth results, with some (~7%) asymmetry in 
the post-shock region. Using the spiderweb mesh produces perfectly symmetric solutions except near the center where the deviation from 
symmetry is as large as ~25%. Using the butterfly mesh produces similar mixed accuracy in symmetry. 
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Fig. 17. — Saltzman piston problem. This problem tests the code's ability to resolve shocks that are oblique to the orientation of the 
grid. The top panel shows the initial grid with 100 X 10 zones. The left wall is a piston moving at a constant velocity, 1.0 cm s~^, to the 
right. Initially, the density and internal energy are set equal to 1.0 g cm~^ and 0.0 ergs, respectively, and we use a gamma-law EOS with 
7 = 5/3. The grid and the density colormap are shown in the bottom panel at t = 0.925 s. See ^10.5l for a discussion of the results. 
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Fig. 18. — The Dukowicz piston problem, another test using an oblique mesh. The initial setup involves two regions having an interface 
with a 60° orientation with the vertical. The top panel shows a low-resolution example of the grid used. Region 1 has a density of 1 g 
cm~^ and is resolved with 144 X 120 zones, and region 2 has a density of 1.5 g cm""^ and is gridded with a 160 X 120 mesh, with the vertical 
mesh lines uniformly slanted at 60°. Initially both regions are in equilibrium with P = 1.0 erg cm~'^. The left boundary is a piston with 
a velocity in the positive x direction and a magnitude of 1.48 cm s~^. The piston-driven shock travels from left to right, and encounters 
the interface. The incident shock continues to the lower density region, a transmitted/refracted shock propagates into the higher density 
region, a vortex sheet develops behind the transmitted shoc k, an d a reflected shock propagates into the incident shock's post-shock flow. 
The orientations of the flow are reproduced accurately (see i|10.6l and Table |3}. 
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Fig. 19. — Nonlinear phase of a single mode of the Rayleigli- Taylor instability at t = 12.7.5 s. Tlie top and bottom densities are p\ = 2.0 
g cm~^ and p2 = 1.0 g cm""^, respectively, and the gravitational acceleration points downward with magnitude g = 0.1 cm s~^. The top 
and bottom boundaries are reflecting while the left and right boundaries are p eriodic. From left to right, grid sizes arc 50 X 150, 74 X 222, 
100 X 300 (this is the resolution used in Fig. I10| l. and 150 X 450. See t|10.7l for discussions comparing these results with those in other 
works. 
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Fig. 20. — Single-mode Rayle igh- Taylor instability: investigat ing the artificial viscosity parameter. The same resolution, setup, and time 
are shown as presented in Fig. 1101 and the third panel of Fig. 1191 Each panel shows results for different values of the artificial viscosity 
parameter, ci, that is the coefficient of Cs(V ■ v) (see J7]l. The artificial viscosity parameters are ci = 0.01 (left panel), ci =0.1 (center 
panel) and, ci = 1.0 (right panel). From this analysis, it would seem that one should choose the lower values of ci to accurately represent 
low Mach number Rayleigh- Taylor flows. However, higher values help to eliminate unwanted ringing in post shock flows. 
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Fig. 21. — The interface amplitude vs. time for the single-mode Rayleigh- Taylor instability test. Resolutions and setup are similar to 
those presented in Fig. I2UI We compare the analytic exponential growth rate (solid line) with simulation results (dashed lines) for viscosity 
parameters, ci, of 0.01, 0.1, and 1.0. Simulations with ci = 0.01 and 0.1 manifest exponential growth for several e-folding times, while the 
run with ci = 1.0 seems to follow the linear phase for only 1 s (~ 1/2 e-folding). 
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Fig. 22. — Kelvin-Helmholtz shear instability. The domain is resolved with 256 X 256 zones. The bottom region has = 1.0, and the 
top region has pt = Pb/Xi where x = 8, and both regions have P = 1.0 erg cm""^. A gamma-law EOS is used with 7 = 5/3. The top panel 
shows the development of Kelvin-Helmholtz rolls at t = 5.5 s for tighoar = \cb and tkh 0.523 s. The bottom plot of Fig. 1221 shows the 
interface amplitude (solid line) vs. time and compares to the expected exponential growth (dashed line) for a simulation with fshcar = j(^b 
(green) and f shear = ^<^b (blue). There are three distinct phases in the log-linear plot: an early transient phase, a phase in which the slope 
most closely matches the exponential growth rate, and the subsequent nonlinear phase. 
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Fig. 23. — Density vs. radius for the core collapse of a IS-Mq star. The model is the sl5 model of Weaver & Woosley. Density vs. radius 
for ID (lines) and 2D (crosses) at times 0, 70, 110, 130, 140, and 150 ms after the start of the calculation. Core bounce occurs at 148 ms. 
The 2D grid is composed of a butterfly mesh in the interior with a minimum cell size of ~0.5 km and extends to 50 km. A spherical grid 
extends the domain out to 4000 km. There are 23,750 cells in total, with an effective resolution of ~250 radial and ~100 angular zones. 
The ID grid has 250 zones, with similar resolution to the 2D run at all radii. For the most part, the ID and 2D calculations track one 
another quite well. There is roughly a ~10% difference in the shock radii at 150 ms. 
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Fig. 24. — Radial component of the gravitational acceleration during the core collapse of a 15-Mq star. Similar to Fig. 1231 Note the 
good match between the ID and 2D results. 
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Fig. 25. — Angular velocity evolution during the core collapse of a rotating 15-Mq star. This test problem was calculated on the same 
grid as in Fig. 1231 The initial angular velocity profile is constant on cylinders and is given by n(r) = l/(f + (r/A)^), where A = fOOO km 

and Slo = 2 radians s~^. As expected, the central angular velocity, Qc is proportional to Pc''^ , where pc is the central density. The central 
density compresses from ■~10^'^ g cm~^ at t = ms to ~2.2 xlO" g cm-3 at t = 160 ms. Therefore, at i = 160 ms, Qc should be ~1600 
radians , consistent with results shown in this figure. Other than a slight, but noticeable, glitch at the location of the shock (~180 km), 
the angular velocity evolves smoothly with no evidence of axis effects. 



